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Abstract 


Molybdenum disulfide, MoS, has recently gained considerable attention as a layered material where neighboring layers are only 
weakly interacting and can easily slide against each other. Therefore, mechanical exfoliation allows the fabrication of single and 
multi-layers and opens the possibility to generate atomically thin crystals with outstanding properties. In contrast to graphene, 
C] it has an optical gap of 1.9 eV. This makes it a prominent candidate for transistor and opto-electronic applications. Single-layer 
MOS) exhibits remarkably different physical properties compared to bulk MoS, due to the absence of interlayer hybridization. For 
= instance, while the band gap of bulk and multi-layer MoS, is indirect, it becomes direct with decreasing number of layers. In this 
= review, we analyze from a theoretical point of view the electronic, optical, and vibrational properties of single-layer, few-layer 
and bulk MoS». In particular, we focus on the effects of spinorbit interaction, number of layers, and applied tensile strain on 
the vibrational and optical properties. We examine the results obtained by different methodologies, mainly ab initio approaches. 
T We also discuss which approximations are suitable for MoS% and layered materials. The effect of external strain on the band 
gap of single-layer MoS» and the crossover from indirect to direct band gap is investigated. We analyze the excitonic effects on 
the absorption spectra. The main features, such as the double peak at the absorption threshold and the high-energy exciton are 
presented. Furthermore, we report on the phonon dispersion relations of single-layer, few-layer and bulk MoS». Based on the 
latter, we explain the behavior of the Raman-active Aj, and E), modes as a function of the number of layers. Finally, we compare 
theoretical and experimental results of Raman, photoluminescence, and optical-absorption spectroscopy. 
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1. Introduction 


For many layered materials, it has been established that the 
few-layer or mono-layer phases have distinct properties with re- 
spect to their bulk counterparts. Often these properties are even 
changing between the mono-, bi- and, tri-layer phases. Within 
the layers, the atoms are held together by strong covalent bonds 
while the inter-layer bonds are rather weak and mostly due to 
van der Waals interaction. As a consequence, the layers can 
easily be separated by mechanical exfoliation and single, quasi 
two-dimensional (2D), and few-layer systems of various mate- 
rials can easily be produced.[1] Some examples are graphene, 
hexagonal boron nitride (BN), semiconducting transition metal 
dichalcogenides MX, (M = Mo, W, Ta, and X = S, Se, Te), 
the superconducting metal NbSez, or the elemental 2D systems 
silicene, germanene, and phosphorene [3]. 

Many of these materials have potential for novel technologi- 
cal functionalities. Graphene is the most prominent single-layer 
material [4]. It does not only have outstanding physical proper- 
ties such as high conductivity, flexibility, and hardness (5], but 
it is also a benchmark for fundamental physics. E.g., it displays 
an anomalous half-integer Quantum Hall effect due to the quasi- 
relativistic behavior (linear crossing in the band-structure) of 
the z-electrons[6] [7]. The fascinating properties of graphene 
have paved the way for intense investigations of alternative lay- 
ered materials. 

Electronics and optical applications often require materials 
with a sizeable band gap. For instance, the channel material 
in field-effect transistors must have a sufficient band gap to 
achieve high on/off ratios [8]. In this respect, the semiconduct- 
ing transition metal dichalcogenides (TMDs) can complement 
or substitute the zero-band gap material graphene[9]. Single- 
layer MoS» is an appealing alternative for opto-electronic 
applications with an optical gap of 1.8-1.9 eV, high quan- 
tum efficiency[10|] [11], an acceptable value for the electron 
mobility(12] [13], and a low power of dissipation[14]. It has 
potential application in nanoscale transitors [9] [15] [12] [8], pho- 
todetectors [16][17], and photovoltaics applications (18][19]|20]. 
Other TMDs such as single-layer WS» also exhibit high photo- 
luminescence yield [21]. 

In this stimulating scenario, TMDs are being intensively in- 
vestigated. Fabrication techniques such as the mechanical ex- 
foliation and the liquid exfoliation produce single- 
and multi-layer crystals with high crystalline quality at low 
cost. This has increased notably the amount of research groups 
working in both fundamental and applied aspects of TMDs. 
Concerning the electrical and optical properties of single-layer, 
multi-layer and bulk MoS», extensive experimental investiga- 
tions have been carried out within the last few years. The most 
important techniques are photoluminescence, optical absorp- 
tion, and electroluminescence spectroscopy [10] [11] 25] 26]. It 
is widely accepted that single-layer MoS» has a direct band gap 
that transforms into an indirect gap with increasing number of 
layers. Similarly, bandgap engineering is possible by applying 
strain. The application of strain drives a direct-to-indirect band 
gap transition in single-layer MoS» [32]. 


Moreover, suitable hydrostatic pressure reduces the band gap 


of single- and multi-layer MoS, resulting in a phase transi- 
tion from semiconductor to metal (33][34]. The group symme- 
try and the spin-orbit interaction in MoS) also raise interesting 
properties. The control of the valley polarization of the photo- 
generated electron-hole pairs paves the way for using MoS» in 
applications related to next-generation spin- and valleytronics 
(35) (36) 37] [B8] [39]. Further studies dealing with charged exci- 
ton complexes (trions) [40] 41] or with second harmonic gener- 
ation have also been published [42] [43]. 

Many challenges remain to be solved in the field of TMDs. 
The problem of obtaining high hole mobility in single-layer 
MoS, hinders the realization of p-n diodes. A proposed so- 
lution is using a monolayer WSez diode, in which the p-n junc- 
tion is created electrostatically by means of two independent 
gate voltages [45]. Another active research field is the 
design of Van der Waals heterostructures. Assembling atomi- 
cally thin layers of distinct 2D materials allows to enrich the 
physical properties [46]. Techniques like chemical vapor depo- 
sition and wet chemical approaches are triggering the fabrica- 
tion of heterostructures 48]. For example, flexible photo- 
voltaic devices of TMDs/graphene layers exhibit quantum ef- 
ficiency above 3096 [49]. A photovoltaic effect has also been 
achieved using a MoS?/WSe» p-n heterojunction [19]. The dif- 
ferent stacking configurations and the band alignments are im- 
portant aspects in bilayer heterostructures[50] [51] [52]. 

Another important activity in the field of TMDs is the charac- 
terization of the vibrational properties of MX». Earlier studies 
of bulk MoS, using Raman and infrared spectroscopy 
as well as Neutron scattering[55] and electron-energy-loss 
spectroscopy[56] had already well characterized the phonons 
at I and the phonon dispersion. In the recent years, a large 
number of Raman studies on mono- and few-layer systems has 
emerged [65]. The Raman fre- 
quencies are correlated with the number of layers which al- 
lows their unequivocal identification. The trend of the Raman 
modes E», (in-plane mode) and Aj, (out-of-plane mode) with 
the number of layers has been intensively discussed, both theo- 
retically and experimentally [68]. The 
Aj, mode follows a predictable behavior. Its frequency grows 
with increasing number of layers, due to the interlayer interac- 
tion. The E), mode shows the opposite trend, i. e., decreasing 
in frequency for an increasing number of layers. 

The experimental findings are accompanied by a vast theo- 
retical literature. The characteristic stacking of ultra-thin layers 
of MoS. adds new challenges to the theoretical approaches. 
For example, the layer thickness influences the dielectric con- 
stant which becomes strongly anisotropic. This enhances the 
Coulomb interaction between carriers. The calculation of the 
excitations has to include these dimensional effect for apply- 
ing accurately the GW method and the Bethe-Salpeter equation. 
For instance, one consequence is an exciton binding energy in 
some layered materials of hundreds of meV, much larger than 
in bulk semiconductors. Also, new models have been devel- 
oped to explain the interplay of the spin-orbit interaction and 
the crystal symmetry (which is layer dependent), and its con- 
sequences, like the valley-Hall effect. Moreover, the interlayer 
interaction has demanded the improvement of the modelling of 


the van der Waals interaction in extended systems. The pre- 
cision of the Raman spectroscopy has allowed to evaluate the 
accuracy of ab initio methods for calculating phonons, and it 
proves how useful is the simulation of the vibrational proper- 
ties for understanding the interlayer interaction and the chem- 
ical bonding. Therefore, the research on layered materials has 
contributed to the appeareance of new methods and to the refor- 
mulation of existing ones. In this review, we give an overview 
of the challenges in the modelling of the spectroscopic prop- 
erties of MoS» and the solutions proposed. The discussion 
of the literature results is complemented by additional calcula- 
tions. We believe the topics discussed here will be also useful 
in the modelling and understanding of other two-dimensional 
materials. 


(a) 


Single-layer 


(b) 


Figure 1: (Color online) (a) MoSzbulk and single-layer. The interlayer distance 
is denoted by d (distance between Mo atoms of different layers). (b) Top view 
of the MoS? single-layer unit cell. 


2. Structural properties 


Bulk molybdenum disulfide (MoS;) belongs to the class of 
transition metal dichalcogenides (TMDs) that crystallize in the 
characteristic 2H polytype. The corresponding Bravais lattice 
is hexagonal and the space group of the crystal is P63/mmc 
(Dg, non-symmorphic group). The unit cell is characterized by 
the lattice parameters a (in-plane lattice constant) and c (out-of- 
plane lattice constant). The basis vectors are 


(la, - a, 0), 


aj = 
a, = (la, ¥a,0), (1) 
a= (0, 0,c). 


The unit cell contains 6 atoms, two Mo atoms are located at the 
Wyckoff 2c sites and four S atoms at the Wyckoff 4f sites. With 
the internal parameter z, the positions, expressed in fractional 
coordinates, are +(1/3, 2/3,1/4) for the Mo atoms and +(2/3, 
1/3,z) as well as +(2/3, 1/3,1/2-z) for the S atoms. 

The single-layer contains one Mo and two S atoms. In this 
case, the inversion symmetry is broken and the space group 
(more precisely, layer group) is P6m2 (Di , Symmorphic group). 
The double-layer is constructed by adding another S-Mo-S 
layer, having now the layer group P3ml (D3 4 Symmorphic 
group). Consequently, an odd number of layers has the same 
symmetry as the single-layer (absence of inversion symmetry), 
whereas an even number has the symmetry of a double-layer 
(with inversion symmetry). 

The crystal structure of MoS»can be specified as a stacking 
of quasi-two-dimensional (2D) S-Mo-S layers along the c di- 
rection. Within each layer, Mo atoms are surrounded by 6 S 
atoms in a trigonal prismatic geometry as illustrated in Fig. 
The bonding type is predominantly covalent within the atom- 
ically thin S-Mo-S layers, whereas the layers themselves are 
weakly bound by Van der Waals (VdW) forces in the crystal. 
The inherent weakness of the interlayer interactions can result 
in different stacking sequences and therefore in different poly- 
typisms as shown in Ref. [74]. 


Defining the optimized geometry is the first step for any cal- 
culation of the phonon spectra and/or the band structure. Most 
of the previous investigations used density-functional theory 
(DFT) on the level of the local-density approximation (LDA) or 
the generalized-gradient approximation (GGA) [66]. We want 
to emphasize that in DFT the accuracy of the calculated quanti- 
ties is determined by the treatment of the exchange correlation 
(XC) energy given by the XC functional. However, the stan- 
dard local (LDA) and semilocal (GGA) XC functionals do not 
account for the long-range van der Waals interactions, which 
are responsible for the stable stacking of the layers and thus 
particularly relevant in two-dimensional materials. Neverthe- 
less, the well-known LDA overestimates the (weak) covalent 
part of the interlayer bonding and compensates thus the miss- 
ing vdW forces yielding a bound ground state for most layered 
materials. This explains the success of LDA in obtaining the ge- 
ometry of many layered materials such as graphite [75], boron 
nitride or graphene on different substrates [78][79] [80]. 
The good performance of LDA in layered materials (although 
fortuitous) has made this approximation widely used in the cal- 
culations of structural properties. 

In the present work, the calculations were partly performed 
with the Vienna ab initio simulation package (VASP)(8T] 
utilizing the projector augmented plane wave (PAW) method 
to describe the core-valence interaction. PAW poten- 
tials with non-local projectors for the molybdenum (Mo) 45, 
4p, 4d, 5s as well as sulfur (S) 35, and 3p valence states were 
generated to minimize errors arising from the frozen core ap- 
proximation. The valence electrons were treated by a scalar- 
relativistic Hamiltonian and spin orbit coupling (SOC) was 
self-consistently included in all VASP calculations as described 
elsewhere [85]. VASP uses DFT with a variety of XC func- 
tionals ranging from LDA to different types of GGAs, to hybrid 
functionals, and VdW density functionals. Furthermore, VASP 
has an implementation of many body perturbation theory such 
as the GW approximation ranging from the single-shot Go Wo 
to a selfconsistent GW (scGW) approximation[87] [88]. 
Concerning standard DFT results presented in this work, the 
XC energy was treated within the LDA(89] and the GGA. For 
the latter, the parametrization of Perdew, Burke, and Ernzerhof 
(PBE), in particular the PBEsol functional was used. 

In order to improve the theoretical lattice parameters cal- 
culated within DFT-LDA/GGA[66], we have also studied the 
structural properties including VdW interactions starting from 
the experimentally observed structural parameters summarized 
in Tab.|1| For this purpose we used the optB86b-VdW func- 
tional, recently implemented in VASP [91]. The optB86b func- 
tional is a non-local correlation functional that approximately 
accounts for dispersion interactions. It is based on the VdW-DF 
proposed by Dion et al. [92], but employs an accurate exchange 
functional particularly optimized for the correlation part. 

The structural optimization of hexagonal bulk MoS» requires 
a 2D energy minimization, since the ground state energy de- 
pends on two degrees of freedom, i. e., the volume and the 
c/a ratio. The experimentally observed structural parameters 
summarized in Tab. [I|have been used as starting point for the 
calculation of the electronic properties of bulk as well as sin- 


Table 1: Bulk MoS? experimental lattice parameters a, c, internal parameter z, and bulk modulus B. 


a (A) c (À) z cla B(GPa) 
Ref. 3.160 12.294 0.621 3.890 
Ref. 3.161 12.295 0.627(5) 3.890 
Ref. 3.140 12.327 3.926 53.4+1.0 
Ref. 3.168(1) 12.3221) 0.625 3.890 
gle layer (1L) MoS, within DFT and methods that go beyond 12.8 


as described in detail below. This minimization was performed 
manually using LDA and PBEsol as well as the optB86b-VdW 
functional, for comparison by varying the unit cell volume of 
bulk MoS. within +10% of the experimentally observed equi- 
librium volume (Vo). For each chosen volume the c/a ratio 
was first optimized by fitting the energy dependence on c/a to 
the Murnaghan equation of state (EOS). The final DFT- 
optimized unit cell volume was obtained by subsequently fitting 
E(V) to the Murnaghan EOS. Note that in each single optimiza- 
tion step the atomic positions were also relaxed by minimizing 
the total forces on the atoms until they were converged to 0.05 
eV/A. 

In order to avoid effects from the changes in size of the ba- 
sis set due to changes in the unit cell volume V, the kinetic en- 
ergy cutoff Ecu has been increased to 350 eV. Convergence with 
respect to k sampling within the Brillouin zone was reached 
with 16 x 16 x 16 I-centered meshes in case of optB86b- 
VdW and with 12 x 12 x 12 I-centered meshes for LDA and 
PBEsol. The manually performed structural optimization was 
cross checked with VASP calculations employing minimization 
algorithms parallel for the atomic positions and the c/a ratio for 
selected volumes in the range +5% of Vo and one subsequent 
Murnaghan EOS fit. From these calculations the Bulk modulus 
is obtained from " 

B=- Lge : (2) 
V av? V-V 

In Table [2| the results of the structure optimization corre- 
sponding to different functionals are summarized. When using 
LDA or optB86b-VdW functionals, the theoretical values of B 
for bulk MoS, agree well with the experimental value given 
in Table [I] Concerning lattice parameters, we observe a small 
underestimation of the in-plane parameter a, both in LDA (1.3 
%) and PBEsol (0.7 %). Contrary, the c parameter is underesti- 
mated by 1.6% and overestimated by 2.4% in LDA and PBEsol, 
respectively. Compared to LDA, the PBEsol functional yields 
the larger deviation of the resulting c/a ratio. The improve- 
ment after including the VdW interactions is substantial, with 
an error of only 0.09 % (see Fig. B]for a comparison with the 
experimental results). 

In conclusion, van der Waals functionals give the most accu- 
rate results for lattice parameters and the bulk modulus. LDA 
tends to underestimate the interlayer distance and the c parame- 
ter, but in average gives acceptable results and it can be trusted 
in the prediction of structural properties. 

Based on the ground state structures summarized in the bulk 
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Figure 2: Comparison of the theoretical a and c lattice parameters, obtained by 
LDA, PBEsol and vdW-DF, with the experimental values of Refs. 


MOS) charge density was calculated on a I'-centered 12x12x3 
k mesh by converging the total energy to 0.1 meV using a ki- 
netic cutoff energy of 350 eV and a Gaussian smearing with a 
smearing width of 50 meV. Tests with 18x18x3 k point grids 
have shown that the electronic band gaps are converged within 
20 meV compared to the results obtained with the 12x12x3 
grid. 

The single-layer MoS, structure has been constructed from 
the optimized bulk structure (Tab. |2) by selecting only the bot- 
tom S-Mo-S layer and adding 20 A vacuum along c direction. 
The atomic positions in the slab geometry have again been re- 
laxed (force convergence criterion of 0.05 eV/À) before calcu- 
lating the band structure on I-centered 12x12x1 K-point grids. 
Convergence tests of the eigenvalues as a function of the vac- 
uum space between repeated layers has been performed up to 
an accuracy of 15 meV, establishing a convergence distance of 
20 À. 

As we will demonstrate later, the main features of the band 
structure of MoS; critically depend on the lattice optimization 
and even small differences can induce significant changes. In 
the case of the phonon band structure, deviations are reflected 
in a rigid shift of the phonon frequencies but the main trends are 
less affected than in the case of the electronic band structure. 


Table 2: Structural parameters of bulk MoS» obtained by minimizing E(V,c/a) with different XC functionals. a and c denote the lattice constants, z the internal 
parameter specifying the atomic positions, B the bulk modulus, and d the interlayer distance defined according to Fig. |1| The uncertainty stemming mainly from 


the EOS fitting in a, c, and c/a is +0.001 A, +0.01 A, and +0.01 A, respectively. 


a(À) c(A) z cja B(GPa) d(À) VdW gap (A) 
LDA 3.120 12.00 0.1214 3.87 40-43 6.039 2.933 
PBEsol 3.138 12.60 0.1264 4.01 18-21 6.305 3.188 
optB86b-VdW 3.164 12.40 0.1236 3.92 39-40 6.203 3.068 


3. Lattice dynamics of MoS, 


The knowledge of the phonon dispersion in a material is in- 
dispensable for the understanding of a large number of macro- 
scopic properties such as the heat capacity, thermal conduc- 
tivity, (phonon-limited) electric conductivity, etc. Vibrational 
spectroscopy (Raman spectroscopy and Infrared absorption 
spectroscopy)[54] give access to the phonons at the Brillouin 
zone center (T point). Inelastic neutron scattering citeWak- 
abayashil975 allows to measure (almost) the full phonon dis- 
persion. Precise semi-empirical modeling of the phonon dis- 
persion and ab-initio calculations in comparison to experimen- 
tal data are a challenge by itself. However, precise modeling 
is also required because details in the vibrational spectra may 
also carry some information about the number of layers and the 
underlying substrate. For graphene, this has been widely ex- 
plored: the so-called 2D line in the spectra splits into sub-peaks 
when going from the single to the multi-layer case [94] [95]. Last 
but not least, the 2D-line also changes position as a function 
of the underlying substrate[96] [80]. All these features 
are related to the double-resonant nature[99] of Raman scatter- 
ing in graphene and on the dependence of the highest optical 
mode on the screening. For MoS» (and related semiconduct- 
ing transition-metal dichalcogenides), the layer-dependence of 
the vibrational spectra is less spectacular than for graphene. 
Nevertheless a clear trend in the lower frequency inter-layer 
shear and breathing modes with increasing layer number can 
be observed((60] (similar to graphene [100]). Also in 
the high-frequency optical modes at T, a clear trend from single 
to multi-layer MoS, has been observed(57] and reproduced in 
other dichalcogenides ; 

Interestingly, the behaviour of the phonon frequencies as a 
function of the number of layers does not always follow the 
intuitive trend. For instance, the frequency of the Raman ac- 
tive in-plane Ej, phonon decreases with the increment of the 
number of layers. This is contrary to the expectation that the 
weak interlayer forces should increase the restoring forces and 
consequently result in an increase of the frequency. The Ra- 
man active A;, does follow the expected behaviour, increasing 
the frequency with the number of layers. Several attempts have 
been done to explain this trend [671162]. This will be criti- 
cally reviewed in this section. 


General Theory 


Before discussing the phonons of bulk and few-layer MoS», 
we briefly review the ab-initio calculation of phonons. (A com- 
plete discussion of the theory of phonons can be found, e.g., 
in Ref. (102].) Starting from the equilibrium geometry of the 
system (see details in Section B}. one obtains the phonon fre- 
quencies from the solution of the secular equation 


1 E 
= Cia) - e (q)| = 0, (3) 
VMM; 
where q is the phonon wave-vector, and M; and M; are the 
atomic masses of atoms J and J. The dynamical matrix is de- 
fined as 


OE 


es (4) 
du’ (q)àw (q) 


Cis jg (q) = 
where u7(q) denotes the displacement of atom Z in direction 
a. The second derivative of the energy in Eq. [4] corresponds 
to the change of the force acting on atom / in direction œ with 
respect to a displacement of atom J in direction 6 [102]. The 
elements of the dynamical matrix at a given wave-vector q can 
be obtained from an ab-initio total energy calculation with dis- 
placed atoms in a correspondingly chosen supercell (that needs 
to be commensurate with q). Another approach consists in the 
use of density functional perturbation theory (DFPT)[103][104], 
where atomic displacements are taken as a perturbation poten- 
tial and the resulting changes in electron density and energy 
are calculated self-consistently through a system of Kohn-Sham 
like equations. Within this approach the phonon frequency can 
be obtained for any q while using the primitive unit-cell. An- 
other way to obtain the dynamical matrix is to use empirical 
interatomic potentials or force constants. A decent fit of the 
MoS, phonon dispersion from Stillinger- Weber potentials has 
recently been suggested by Jiang et al.[105]. Such a fit allows 
to study much larger systems such as nanoribbons, nanotubes 
or heterostructures. The advantage of DFPT is, however, the 
higher accuracy and the automatic inclusion of mid and long- 
range interactions. For later reference, we note that real-space 
force constants can be obtained from the dynamical matrix by 
discrete Fourier Transform: 


Croan = 5 Y et^ Cu. (5) 
where N, is the number of unit-cells (related to the density of 
the q-point sampling). The physical meaning of the C; ;g(R) 
is the force acting on atom / in direction a as atom J in a unit 
cell at R is displaced in direction £ and all other atoms in the 
crystal are kept constant. 

It is important to note that for polar systems, the phonons in 
the long-wavelength limit can couple to macroscopic electric 
fields. In mathematical terms, this means that the dynamical 
matrix in the limit (q — 0) can be written as the sum of the dy- 
namical matrix at zero external field and a "non-analytic" part 
that takes into account the coupling to the electric field and de- 
pends on the direction in which the limit q — 0 is taken: 


Cj, jg(q > 0) = Čia ;g(q = 0) + Cr (Qd 2 0). (6) 


The non-analytic part contains the effect of the long-range 
Coulomb forces and is responsible for the splitting of some 
of the longitudinal optical (LO) and transverse optical (TO) 
modes. [103] [104]: Its general form is as follows: 


An 2 (q g Z* 1)e(q g Z* yg 


CYA B = Q q-e-q 


; (7) 
where Q is the volume of the unit cell, Z*, stands for the Born 
effective charge tensor of atom s and e is the dielectric tensor. 
Since the dielectric tensor is fairly large in bulk MoS% (e = 


6y = 15.4, & = 7.43, the effect of LO-TO splitting is visible, 
but not very pronounced (< 2.6 cm 1). 

We will discuss in the following first the phonon dispersions 
of bulk and single-layer MoS». Afterwards, we will discuss in 
detail the symmetry of bulk and single-layer phonons at I. We 
will single our the Raman and infrared (IR) active modes and 
discuss how their frequencies evolve with the number of layers. 


Phonon dispersion 


The phonon dispersions of single-layer and bulk MoS% are 
shown in Figure [3] Overall, the single-layer and bulk phonon 
dispersions have a remarkable resemblance. In the bulk, all 
single-layer modes are split into two branches but since the 
inter-layer interaction is weak, the splitting is very low (similar 
to the situation in graphite and graphene.[75] The only notable 
exception from this is the splitting of the acoustic modes of bulk 
MoS) around T. 

We have also depicted the experimental data obtained with 
neutron inelastic scattering spectroscopy for bulk MoS»[55] as 
well as the result of IR absorption and Raman scattering at T. 
The overall agreement between theory and experiment is rather 
good, even for the inter-layer modes. This confirms our expec- 
tation that the LDA describes reasonably well the inter-layer 
interaction (even though not describing the proper physics of 
the inter-layer forces, as discussed in Section[2). 

The bulk phonon dispersion has three acoustic modes. Those 
that vibrate in-plane (longitudinal acoustic, LA, and transverse 
acoustic, TA) have a linear dispersion and higher energy than 
the out-of-plane acoustic (ZA) mode. The latter displays a q?- 
dependence analogously to that of the ZA mode in graphene 
(which is a consequence of the point-group symmetry [108]). 
The low frequency optical modes are found at 35.2 and 57.7 
cm! and correspond to rigid-layer shear mode and layer- 
breathing mode (LBM) respectively (see left panels of Fig.|5) in 
analogy with the low frequency optical modes in graphite[75]). 
It is worth to mention the absence of degeneracies at the high 
symmetry points M and K and the two crossings of the LA 
and TA branches just before and after the M point. The high 
frequency optical modes are separated from the low frequency 
modes by a gap of 49 cm !. 

The single-layer phonon dispersion is very similar to the bulk 
one. The number of phonon branches is reduced to nine. At low 
frequencies, the shear mode and layer-breathing mode are ab- 
sent. At higher energies, very little difference between bulk and 
single-layer dispersion can be seen. This is due to the fact that 
the inter-layer interaction is very weak. The subtle splitting and 
frequency shifts of zoner-center modes in gerade and ungerade 
modes (as going from single layer to the bulk) will be discussed 
below. 

The densities of states (DOS) of single-layer and bulk are 
represented in the right panels of Fig. The differences be- 
tween single-layer and bulk DOS are minimal, except a little 
shoulder around 60 cm ! in the bulk DOS due to the low fre- 
quency optical modes. In both cases the highest peaks are lo- 
cated close to the Raman active modes E), and Aj, and they are 
due to the flatness of the bands around I. 


Figure 4: Raman spectrum from Ref. [58], recorded using 623.8 nm excitation 
(red dots). One-phonon density of states (blue lines) and two-phonon density 
of states (green lines). 


Intensity (arb. units) 


The density of states can be partially measured in 2nd and 
higher-order Raman spectra. We have represented in Fig. 
the Raman spectrum of MoS» bulk of Ref. (58], obtained by 
exciting the sample with a laser frequency of 632 nm. Similar 
spectra can be found in older studies and other re- 
cent studies [106][107]. First order Raman peaks are due to the 
excitation of zero-momentum phonons. We can identify in the 
spectrum the modes E), and A;,. Moreover, working in condi- 
tions of resonant Raman scattering, second-order Raman modes 
can be obtained. These modes come from the addition or sub- 
traction of modes with opposite momentum, w;(q) + w;(q), to- 
gether with the resonance of an intermediate excited electronic 
state[101]. The result is a rich combination of Raman modes, 
as shown in Fig. |3| The identification can be done with the help 
of the density of states. We have calculated the 1-phonon and 
the 2-phonon density of states for MoS» bulk, represented by 
the solid blue line and the dashed green line, respectively. 

A careful examination of the 2-phonon density of states tells 
us which phonons are participating in the Raman spectra. Thus, 
the longitudinal acoustic branch at M, denoted as LA(M), cou- 
ples to the modes Ej, E» and A;,, resulting in overtones in 
the Raman spectrum. Other combinations include 2 x E), or 
2 x Ajg, always with momentum M. The second-order Raman 
modes are much more restrictive than first-order modes. The 
concurrence of phonon modes and electronic levels is needed. 
Such alignment depends strongly on the electronic structure. 
Consequently, the second-order Raman spectrum has revealed 
useful to establish the fingerprints of single-layer systems with 
respect to the bulk, as discussed in Ref. [101]. 


Symmetry analysis of phonon modes 


We have drawn in Figs. |5| and [6] the atomic displacements 
(eigenvectors) of optical phonon modes of bulk and single-layer 
MoS, at I. Group theoretical analysis[53] 
yields for the 15 optical modes of bulk MoS; (Dg, symme- 
try) the following decomposition in irreducible representations: 
Aig P Ary, ® 2B», ® Biu ® Eig ®@ Eiu ® 2Ezg @ E>,. The Aig, Eig 
and E», modes are Raman active and the A2, and £,, modes are 
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Figure 3: Phonon dispersion curves and density of states of (a) single-layer and (b) bulk MoS2. The symbols denote experimental data. Blue circles: neutron 
scattering , yellow circles: first and higher-order Raman scattering of bulk MoS2 , red dots Raman scattering of single-layer MoS» : 
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Figure 5: Sketch of the optical phonon modes of bulk MoS. In the first row, the modes with polarization (atom-movement) parallel to the layers are plotted in 
ascending order. In the second row, the perpendicular modes are shown. “Davydov pairs" of phonon modes are plotted in one box. The phonon frequencies (in 
cm7!) are the calculated values of Ref. [67]. 
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Figure 6: Sketch of the optical phonon modes of single-layer MoS?. The phonon frequencies (in cm7!) are the calculated values of Ref. [67]. 


infrared active. For the 15 optical modes of double-layer MoS» 
(D34 symmetry), the decomposition is: 341, 6245, 63E, 6 2E, 
where the gerade modes are Raman active and the ungerade 
ones are IR active. For the 6 optical modes of the single- 
layer, one obtains the following irreducible representations: 
A, € A5 e E' GE”. The Aj and E" modes are Raman active, 
the E’ mode is both IR and Raman active. 


The attribution of the different symmetries to the phonon 
modes in Figs. [5]and[6|can be understood quite intuitively (tak- 
ing into account that the drawings are simplified 2D versions 
of the real 3D modes). All E modes are doubly degenerate and 
correspond thus to in-plane vibrations of Mo and/or S atoms be- 
cause a rotation by any angle yields another phonon-mode with 
the same frequency. The non-degenerate A and B modes must 
therefore correspond to perpendicular movement of the atoms. 
For bulk MoS;, (space group P63/mmc, point group Den), there 
is an inversion center half-way between two S atoms of neigh- 
boring layers. We can thus distinguish between gerade (g) and 
ungerade (u) modes which are symmetric/anti-symmetric with 
respect to inversion. The gerade modes are those where atoms 
in neighboring layers move with a phase shift of z while the 
ungerade modes correspond to in-phase movement. All phonon 
modes of bulk MoS» thus come in pairs of gerade and unger- 
ade modes which are close in frequency. This can be clearly 
seen for the modes in panels 1 to 4 of Fig. Furthermore, 
the “shear mode" at 35.2 cm™! is the gerade counterpart of the 
in-plane acoustic mode (not shown) and the “layer-breathing 
mode" (LBM) at 55.7 cm”! is the gerade counterpart of the out- 
of-plane acoustic mode. In almost all cases, the gerade mode 
is higher in frequency than the ungerade mode. This is because 
the weak (Van-der-Waals like) bond between S atoms of neigh- 
boring sheets is elongated and squeezed in the gerade mode 
(thus gives rise to an additional restoring force) but kept con- 
stant in the ungerade mode. The notable exception is the case of 
the modes in panel 2 where the ungerade mode is higher in en- 
ergy. We will come back to this important case in the next sub- 
section. One can easily see that only ungerade modes can be IR 
active: for a mode to be IR active, a net dipole must be formed 
through the displacement of positive charges in one direction 
and negative charges in the opposite direction. However, in ger- 
ade modes, the dipoles formed on one layer are canceled out by 
the oppositely oriented dipoles on the neighboring layer. Since 


in systems with inversion symmetry, a phonon mode cannot be 
both IR and Raman active, only the gerade modes can be Ra- 
man active in bulk MoS;. 

The distinction between A and B modes is made by rotat- 
ing the crystal by 27/6 around the principal rotation axis. This 
rotation is a non-symmorphic symmetry, i. e., it has to be ac- 
companied by a translation normal to the layer-plane in order 
to map the crystal into itself. In our reduced 2D representation 
of the vibrational modes this corresponds to a translation of the 
3 atoms of the upper layer onto the 3 atoms of the lower layer. 
If the arrows change direction, the mode is B, otherwise A. Fi- 
nally, for the singly degenerate modes, the subscript 1 (2) stand 
for modes that are symmetric (antisymmetric) with respect to 
rotation around a C» axis crossing an Mo atom perpendicularly 
to the 2D plane of projection. For the doubly degenerate E 
modes, it is the other way around. 

For even N-layers of MoS», the space-group symmetry is 
P3m] and the assignment of the phonon-mode symmetries has 
to be done according to the D34 point-group symmetry. Since 
inversion symmetry is present, the mode assignment is very 
similar to the one of bulk MoS». For the doubly degenerate 
E modes (see Fig. 5), the subscripts 1 and 2 are dropped. All 
E,, modes are IR active and all E, modes are Raman active. Out 
of the perpendicularly polarized modes, the inactive Bj, mode 
turns into an IR active A>, mode, the inactive B) modes turns 
into a Raman active Aj, modes. Notably, the layered breathing 
mode (LBM) is, in principle, Raman active. Indeed, for dou- 
ble and 4-layer MoS, this mode has been detected in Raman 
measurements, albeit with small amplitude|[64]. 

For the single layer and for add-numbered multi-layers, the 
space group is P6m2 and the corresponding point-symmetry 
group is D3,. Since inversion symmetry is absent in this group, 
there is no distinction between gerade and ungerade modes. In- 
stead, modes that are symmetric under c; (reflection at the xy- 
plane) are labeled with a prime and anti-symmetric modes with 
a double prime (Fig. (6). 

The experimental and theoretical frequencies of all phonon 
modes of single-layer and bulk MoS» at I are summarized in 
Table [3] For the IR active modes of bulk MoS», we give both 
the values for longitudinal-optical (LO) and transverse-optical 
(TO) modes. The LO-TO splitting is calculated from the non- 
analytic part of the dynamical matrix (Eq.[7) which only affects 


Table 3: Phonon modes at T of bulk and single-layer MoS» (inspired by Table II of Ref. (54]). The polarization of the modes is in-plane (xy) or perpendicular 
z). The irreducible representation (Irrep.) of each mode is calculated from the corresponding point-symmetry group (Dg; for bulk, D3; for single-layer). For the 
character of the modes, we distinguish between Raman active (R), infrared active (IR), acoustic modes (a), and inactive modes (i). 


single layer (D3, symmetry) bulk (Dg, symmetry) 
Pol. Atoms | Irrep. Char. Cale £3 M eg Irrep. Char. cute f] RM 
xy Ms | E a 0.0 E E B 33 
z Mos | A; a 0.0 Br : Es 
xy s |E' R 289.2 E n oe 287 
E, R 387.8 383 
xy Mos | F’ R+IR(E Lc) 391.7 384.3 E, IRE Lc) T n vie 
Z S Ai R 410.3 403.1 ij : : R 409 
z Mos | A; IR(Elle) 476.0 Am TRUM) i "e 3 p 
B5, i 473.2 


the IR active modes. 


Anomalous Davydov splitting 


As mentioned above, in bulk MoS;, all modes appear as pairs 
of gerade and ungerade modes (Fig.[5]and Table B). This is be- 
cause the unit cell of bulk MoS, contains 6 atoms while the 
single-layer unit cell only contains 3. The frequency splitting 
of gerade and ungerade modes in layered materials and molec- 
ular crystals is known as “Davydov splitting" or “factor group 
splitting” [113]. The model of Davydov has been used 
in particular to explain the splitting of modes that are IR and 
R active for the single-layer into a Raman active mode and an 
IR active mode of the bulk (mode No. 2 in Figs. 5]and[6). It 
was recognized early that neither the weak Van-der- 
Waals coupling between neighboring layers nor a simple model 
of dipolar couplings matches the experimental observation that 
YRaman < VIR,TO < YIR,LO for some layered materials and, 
in particular, MoS. A model involving quadrupole interaction 
was proposed by Ghosh et al.[[115][1 16] but could not be under- 
pinned by numerical calculations. 

The explanation of the “normal” Davydov splitting in van- 
der-Waals bonded layered materials is straightforward. As can 
be seen in Fig.|5| the weak inter-layer bonding can be viewed 
as an additional (weak) spring constant acting between sul- 
fur atoms from neighboring layers (red dashed lines). For the 
ungerade modes, the S-atoms are moving in phase and the addi- 
tional spring thus is not “used”. However, for the gerade modes, 
where the S-atoms are vibrating with a phase shift of x, the ad- 
ditional spring is elongated and compressed and thus yields an 
additional restoring force. This leads, in general, to an upshift 
of the frequencies of the gerade modes. Since the interaction 
is weak, the frequency shift is small (< 5cm^!). Furthermore, 
the effect is more pronounced for the perpendicularly polarized 


modes than for the in-plane modes (Fig. [5]and Table 3). The 
only exception from the “normal” Davydov splitting is thus the 
mode No. 2. One might argue that this case is exceptional, 
because the LO-shift of the E} mode makes its frequency 
higher then the one of the Ej, mode. However, even without 
the LO-shift, experiments[54] and calculations[67] agree that 
YE.) + Y(E,, TO). 


Since our ab-initio calculations reproduce the anomalous 
Davydov splitting for mode No. 2, we can use the interatomic 
force constants, derived from the calculations, in order to find 
the origin of this seemingly anomalous behaviour. The situa- 
tion is depicted in Fig.[7](a) and (b). We analyze in detail the 
force constants between S atoms of neighboring layers (blue 
springs) and also the force constants between S atoms on one 
layer with Mo atoms on the neighboring layer (red springs). In 
the Ej mode it is the sum of all the S—S spring constants that 
leads to an additional restoring force and thus to an up-shift 
(with respect to the same mode in the fictitious isolated layer). 
However, for the Eiu mode, an additional restoring force arises 
as well, this time due to the Mo-S interactions. As it turns 
out, their effect is stronger than the ones of the S-S springs. 
This follows from the numerical values of the horizontal com- 
ponents of the S-S and Mo-S force constants. We present the 
values in Fig. F]. For a given S atom, we calculate the sum 
of the (horizontal) force constants over all nearest, next-nearest, 
..., S and Mo atoms of the adjacent layer. Negative sign implies 
restoring force (the S atoms is pushed back to the left when 
displaced to the right). In Fig. Fo) one can see that the in- 
teraction of S with the three next-nearest S atoms of the adja- 
cent layer is stronger than the interaction of S with the closest 
Mo atom. However, the cumulative effect of the S—Mo inter- 
actions is larger than the one of the S-S interaction. This ex- 
plains why for this mode pair the sign of the Davydov splitting 
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Figure 7: (a) and (b) Sketch of the Eju and E} modes. (c) Sketch of in-plane force constants ks_s and kmo-s . (d) Phonon frequencies of Eiu (red) and E), (blue) 
modes as a function of the ratio Ky, s /ks_s’. 


is negative. The dominance of the inter-layer S—Mo interaction 
over the inter-layer S—S interaction was already invoked in the 
force-constant model of Luo et al.[62] In that model, all the in- 
teraction was ascribed to the closest atom pairs of adjacent lay- 
ers. This renders the semi-empirical model simple and quantita- 
tively successful. In semi-empirical calculations using the code 
GULPÍ(I117][118], we have verified that it is possible to inverse 
the frequencies of the E), and E, modes by changing the rel- 
ative strength of the cross-layer S-S and S-Mo interaction (see 
Fig.[7|(d)). However, the ab-initio calculations demonstrate that 
the physical reality is more complex. The 2nd-nearest neighbor 
interaction between sulfur atoms across layers is even repulsive. 
Thus, the correct balance between S—S and S—Mo interaction is 
only found by summing over all interactions. It turns out that 
the force constants decrease rather quickly with increasing dis- 
tance and the Coulomb contribution (from the effective charges) 
is rather small. 

Note that the situation is different for mode No. 1 in which 
the Mo atoms do not move and the inter-layer S-Mo interaction 
thus plays a negligible role for the Davydov splitting. In this 
case, the Raman active £;, mode is slightly higher in frequency 
than the inactive E», mode, following the intuitive expectation 
and yielding the “normal” sign of the Davydov splitting. 

We will see in the following that also the dependence of the 
frequencies of the E} on the number of layers follows an un- 
expected trend which can be used for the determination of the 
number of layers via Raman spectroscopy. 


Dependence of Raman active modes on number of layers 

Since the beginning of the research on MoS;flakes, the Ra- 
man modes have been used to identify the number of layers[57| 
[64]. The correspondence between frequency and num- 
ber of layers has been done by comparing with other tech- 
niques such as atomic force microscopy or optical contrast. The 
phonon modes used to characterize the thickness are usually the 
high-frequency Raman modes E), and Aj, (see Fig. 5 and the 
breathing and shearing modes at low-frequency [64]- We will 
summarize in the following the results and analyze the physics 
of the frequency trends. 


High-frequency modes 

In the single layer, the high frequency I' modes El ? and £i, 
collapse into the mode £’. (From Fig. elit is evident that with 
increasing inter-layer distance, the modes Ei and E, acquire 
the same frequency.) Interestingly, as measured in Ref. 
and indicated in Table |3|(see also Pigs, Pfand , the bulk E), 
mode is lower in frequency than the single-layer E' mode. This 
contradicts the expectation that the additional inter-layer inter- 
action should increase the frequency from single-layer to bulk. 
Even the bulk E,,, mode (which is higher in frequency than the 
El mode due to the anomalous Davydov splitting) is slightly 
lower than the single layer E^ mode. The same behaviour (that 
the bulk modes are lower in frequency than the single-layer 
mode) can be observed for the in-plane mode No. 1 (£2, and 
Ej, in bulk versus E" in the single layer) and for the out-of- 
plane mode No. 4 (A5, and B! in bulk versus Aj in the sin- 
gle layer). Only the out-of-plane mode Aj, (No. 3) follows the 


expected trend that the inter-layer interaction increases the fre- 
quency with respect to the single-layer mode Aj. 


Figure [8] shows the frequency of the A; and E’ modes as a 
function of layer number. We compare LDA calculations (cir- 
cles comes from Ref. and triangles from Ref.[62]) with 
the experimental data of Refs. and [62]. The out-of-plane 
mode A! increases in frequency with increasing number of lay- 
ers. The explanation lies in the interlayer interaction, caused 
by weak van-der-Waals bonding of the sulphur atoms of neigh- 
bouring layers. The stacking of layers can thus be seen as the 
addition of a spring between the sulfur atoms of neighboring 
layers. Within the picture of connected harmonic oscillators, 
this results in an increase of the frequency. The LDA calcula- 
tions reproduces well this trend, as can be seen in Fig. [8] The 
small disagreement between theory and experiment can be at- 
tributed to the inaccuracy of the interlayer interaction given by 
LDA. 


The in-plane mode £' displays the opposite trend, decreasing 
in frequency by about 4 cm^! from single-layer to bulk (Fig. 8 
lower panel). This is - at first sight - unexpected, because the 
additional “spring” between the sulfur atoms should lead to an 
increased restoring force and thus to a frequency increase as 
in the case of the A! mode. Several attempts have been made 
in the past to explain this anomalous behaviour, ascribing it to 
long-range Coulomb interactions[57]. In our previous previous 
publication(67], we have investigated how the dielectric screen- 
ing in the bulk environment reduces the long-range (Coulomb) 
part of the force constants. However, the long-range part plays 
only a minor role. We have verified this by performing an ab- 
initio phonon calculation of the E" mode of single-layer MoS» 
sandwiched between graphene-layers. If the distance between 
the sulfur atoms and the graphene layer is higher than 6À, there 
is no "chemical" interaction between the different layers and 
the graphene just enhances the dielectric screening of the MoS2 
layer. Since the E' remains unaffected, we conclude that the 
long-range Coulomb effect can be discarded as a possible effect 
for the anomalous frequency trend. 


The solution to the problem has been given by Luo et al.[62] 
and is related to a weakening of the nearest neighbor Mo-S 
force-constant in the bulk environment. To be precise, one has 
to compare the Real space force constants C Li 5 (0) for the force 
and displacement parallel to the layer. (See Eq. (5) for the defi- 
nition of the force constant). This is the dominant force constant 
determining the frequency of the E" mode as becomes immedi- 
ately clear from the displacement pattern in Fig. [el Force con- 
stants between more distant atoms play a very small role. There 
are two reasons why C Mis 5 (0) is larger for the single-layer than 
for bulk. One reason is that in the single-layer, the Mo-S bond 
is slightly shortened. However, even without this change in 
bond-length, the Mo-S force constant is slightly larger in the 
single-layer than in bulk. This can be obtained in an ab-initio 
calculation of the single-layer using the interatomic distances 
from bulk. The results of our calculations are shown in Table B] 
and are very similar to the values of Fig. 3 in Ref. [62]. The 
frequency of the E" mode is proportional to the force-constant: 


wg X Nous ug (0). Even though the differences seem tiny, they 


Table 4: Ab-initio Force constant cho 5 (0) in (Ha/Bohr?) in single layer (SL) 
and bulk MoS;. 


bulk nó 


-0.1102 


bulk geom. 
-0.1119 


opt. geom. 
-0.1127 


explain the calculated and observed frequency differences be- 
tween single-layer and bulk in table[3] The finding that the Mo— 
S force constant is smaller in bulk than in single-layer, even at 
equal bond-length, is related to a (slight) redistribution of the 
charge density when a neighboring layer is present. 

The fact that the A! and the £’ modes move in opposite di- 
rections with increasing number of layers, makes the distance 
between the two corresponding peaks in the Raman spectra 
a reliable measure for the layer number[57] [62]. But this is 
not the only way to detect the layer number in Raman spec- 
troscopy. The low frequency Raman active modes display an 
even stronger dependence as explained in the following. 

The shearing mode (C), denoted in bulk as ES. is the rigid- 
layer displacement in-plane. This mode is Raman active in 
bulk, as indicated in Table[3] The layer-breathing mode (LBM) 
corresponds to vertical rigid-layer vibrations, in the case of 
bulk, where is has B3, symmetry, it is a silent mode. However, 
in the bi-layer case it has Aj, symmetry and is (weakly) visible. 
Several groups have investigated the low frequency behaviour 
of few-layer MoS> [60][64][63]. The frequency trends as a func- 
tion of the layer number can be explained via a simple analyti- 
cal model that was first developed to explain the corresponding 
modes in few-layer graphene[100) [119]. In this model, N layers 
with a mass per unit area u are connected via harmonic springs. 
One distinguishes between force constants (per unit area) for 
displacement perpendicular œ, and parallel a to the layer, re- 
spectively. Mathematically, the model is equivalent to a linear 
chain of N atoms. Assuming a time dependence Assuming a 
time dependence of u,(t) = u? expliwt] for all the N atoms, 
Newton’s equation of motion lead to the secular equation 


" 1 -1 0 jo Ory, 
-1 2 -I 0 
F _a| 0 ct | 
H -1 0 
0 -1 2 -I 
Un 0 0 0 -1 1 Jj 


(8) 
where œ = «œ, for the layer-breathing modes and œ = ay for 
the shear modes. The frequency of the ith phonon mode (i = 


1,...,N) is 
zi 
Qj = 4| — |1 — cos : 
aE 


For i = 1 one obtains the acoustic (translational) mode, i = 2 
yields the lowest non-vanishing frequency of the mode with a 
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Figure 8: Phonon frequencies of Aj, and E), modes as a function of MoS» layer 
thickness. The symbols corresponds to: (red circle) this work (blue triangles) 
LDA calculations of Ref. [62], (black stars) experimental data from Ref. 
and from Ref. (green diamonds). 


nodeless envelope function, and i = N yields the highest fre- 
quency mode where neighboring layers are vibrating with a 
phase shift of (almost) 7. 

Fig. [9]shows the Raman spectra published in Ref. [64]. The 
number of layers ranges from 1 to 19 in the case of odd number 
of layers (ONL) and from 2 to 18 in the case of even num- 
ber of layers (ENL). Evidently, the single-layer MoS, Raman 
spectra has no signature of low-frequency modes. The peaks at 
4.55 cm ! are due to the Brillouin scattering of the LA mode 
from the silicon substrate. One can see that some peaks stiffen 
(dashed lines) for increasing thickness and others are softened 
(dotted lines). Fig. pro shows the shear (C) and breathing 
(LBM) mode as a function of the number of layers. We can see 
also in Fig. Pra) that the full width at half maximum (FWHM) 
behaves in a different way for the C or the LBM modes. In 
the case of the LBM mode (blue dots in in Fig. p] (a,b)), the 
branch with the largest weight is the one with i = 2. Accord- 
ing to Eq. p] the frequency of this branch as a function of layer 
number N goes like 


o BMN) = w BMO) 1 — cos(x/N), (10) 


where wl BM (2) = «2a,/p. This is the functional form of the 
blue diamonds in Fig. p]. For i = N, the LBM increases in 
frequency according to 


wh BM (yy = 4gBMQ5) 1 + cos(/N) (11) 


and approaches, for N — ov, the value of the Be bulk mode at 


55.7 cm7!. However, since the bulk mode is not Raman active, 
the intensity of this mode quickly decreases with increasing N 
and the mode is already almost invisible for N = 4. For in- 
termediate values of 2 < i < N/2, side branches of the LBM 
appear and are clearly visible in the Raman spectra of Fig. P] 
The same analysis can be done for the shear (C) mode. In 
this case, it is the i = N branch that has the dominant weight 
and converges towards the bulk E, shear mode with increasing 


N: 
WEN) = «C 2) V1 + cos(x/N). (12) 


This is the functional form of of the red squares in Fig. [9] Co. 
The frequency ratio of bulk and double layer shear modes is 
Cbulk/ (2) = V2 which is verified by the experimental data[60] 
[64]. Some side-branches for N < i < N are visible in the 
spectra as well, however with lower intensity than the i = N 
branch. 

Due to the strong layer dependence of the frequencies, the 
shear and compression mode are a very sensitive tool for the de- 
termination of layer-thickness by Raman spectroscopy [601[64]. 
The monoatomic chain model is able to explain the main 
physics of these modes. Small deviations from the analytic re- 
sult have been observed and might be due to anharmonic 
effects[61]. The disadvantage of the use of these modes lies 
in the detection of frequencies close to those of the Brillouin 
scattering. 
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Figure 9: Stokes and anti-Stokes Raman spectra for (a) odd and (b) even number of layers (ONL and ENL, respectively). The spectrum of bulk is included in (a) 
and (b). (c) Raman frequency and (d) FWHM, as a function of the number of layers of the breathing mode (LBM) and shear mode (C) (Reprinted with permission 
from Ref. [64]. Copyright (2013) by the American Physical Society). 
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4. Electronic structure 


This section is devoted to the main properties of the elec- 
tronic structure of MoS». Based on first principles calculations, 
the characteristics of the band structures of single-layer, double- 
layer and bulk MoS»are discussed. In particular, we analyze the 
electronic band gaps and show ab initiohow the band structures 
depend on the unit cell parameters and the structural optimiza- 
tion. We explain the reasons for the differences in the results 
obtained through the different computational approaches. 

Historically, TMDs are used in the field of tribology as lubri- 
cants. The attention given to TMDs decades ago has lead to sev- 
eral theoretical studies of the band structure of MoS, in single- 
layer and bulk forms|120] [21]. These studies were comple- 
mented more recently with angle-resolved photoelectron spec- 
troscopy measurements for bulk MoS. accompanied by ab ini- 
tio calculations [122]. 

The current interest in MoS, [9] [123], the availability of 
high-quality single-layer flakes [24], and the improvement of 
experimental results have prompted new theoretical studies in 
the past 5 years. Regarding the electronic structure, the most 
efficient approach with respect to computational cost and accu- 
racy is the use of DFT-LDA/GGA. Due to the potential applica- 
tion of single-layer MoS, in transistors, most calculations are 
focused on the band gap. By using LDA, single-layer MoS% is 
determined as a direct semiconductor with a band gap of 1.78 
eV at the K point of the Brillouin zone [124]. The transition 
from a direct to an indirect gap semiconductor with increas- 
ing number of layers is also confirmed [27] [28]. The extensive 
use of standard DFT in MoS, (and the comparison with more 
advanced methods) has demonstrated that this approach is ade- 
quate for investigating trends. However, an inherent problem of 
DFT-LDA/GGA is the underestimation of the band gap which is 
related to the lack of the derivative discontinuity in (semi)local 
exchange-correlation potentials[125]. 

In a first attempt to improve the electronic band gaps and 
band dispersions at low computational cost, the modified 
Becke-Johnson (MBJ) potential combined with LDA cor- 
relation was applied. The MBJLDA approach was developed 
by F. Tran and P. Blaha and implemented in VASP [128]. 
The MBJ potential is a local approximation to an atomic exact 
exchange potential plus a screening term (screening parameter 
c) and is given by 


VMBI (py = cV BR) + (3c - 2)- Ta , | ind (13) 


In this expression, pọ denotes the electron density, tT, the 
kinetic energy density, and VPR(r) the Becke-Roussel (BR) po- 
tential [129]. In the parameter-free MBJLDA calculation em- 
ployed in this study, the c parameter is chosen to depend lin- 
early on the square root of the average of |Vp|/p over the unit 
cell volume and is self-consistently determined. 

Alternatively, the screened hybrid functional HSE 
132] and the improved HSEsol functional were tested. 
The success of HSE in predicting band gaps with accuracy 


comparable to that of schemes based on many body perturba- 
tion theory (GW methods) but significantly reduced computa- 
tional cost is multiply demonstrated in the work of G. Scuse- 
ria (see the review in Ref. [134]) as well as in indepen- 
dent studies including a variety of materials and properties 
[85]. In general, hybrid functionals are 
constructed by using the DFT correlation energy E. and adding 
an exchange energy F that consists of 25% Hartree-Fock (HF) 
exchange and 75% DFT exchange. Furthermore, in the concept 
of the screened HSE functional the expensive integrals of 
the slowly decaying long-ranged part of the HF exchange are 
avoided by further separating the E, into a short- (SR) and long- 
ranged (LR) term, where the latter is replaced by its DFT coun- 
terpart. An additional parameter u defines the range-separation 
[139]. The HSEsol functional is analogous to the HSE func- 
tional, but is based on the PBEsol functional for all DFT parts 
according to 
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Concerning the electronic properties, the highest level of ac- 
curacy has been achieved by GW calculations. In this work 
the band structures were studied using the single-shot (Go Wo) 
and the self-consistent GW (scGW) approximation. In both 
approaches, the dynamically (frequency dependent) screened 
Coulomb interaction W and the self energy X(7, r’, w) were cal- 
culated using the DFT-LDA wave functions. In the scGW case, 
both, the quasiparticle (QP) energies (one electron energies) 
and the one electron orbitals (wave functions) are updated in G 
and W. Note that the attractive electron-hole interactions were 
not accounted for by vertex corrections in W. Thus the calcu- 
lated band gaps are expected to be slightly overestimated com- 
pared to experiment [88]. In the particular case of MoS, the 
GW method has been used in many flavours, yielding a signifi- 
cant spread of values for the band gap, to be discussed below. 


4.1. Characterization of the band structure of single-layer and 
bulk MoS» 


First of all, we discuss the main features of the electronic 
structure of MoS. Figure shows the band structure for 
single-layer (1L), double-layer (2L) and bulk MoS». The vdW- 
DF optimized lattice constant (Section) is used in all calcula- 
tions. The suitability of using either the experimental or the the- 
oretical lattice parameter in the calculations is still controver- 
sial. For the calculations presented in this review, we have cho- 
sen the latter, which guarantees a strain-free structure. Thereby, 
we will further be able to investigate strain effects on the elec- 
tronic structure. 

Figure[I0|shows the relevant features of the bulk MoS, band 
structure: 


e two distinct valence band edges (VBEs) located at K and 
T, 


E (eV) 


Figure 10: Band structure of MoS» single-layer (left), double-layer (center), and bulk (right) in the LDA (thin dashed lines) and in the Go Wo approximation (thick 
continuous lines). The lattice parameters have been obtained from the structure optimization using the optB86b-VdW functional. 


e three conduction band extrema (CBEs) at T (half way be- 
tween I and K), K, and È (half way between I' and M), 


e the valence band maximum (VBM) located at F and the 
conduction band minimum (CBM) at T, 


e a fundamental electronic band gap of indirect nature that 
is defined by the energy difference T. — I, 


e the splitting of the VBM at K into states K,; and K,2 due 
to interlayer interaction, 


e two-fold spin degeneracy of all states due to inversion 
symmetry, and 


e nearly parabolic band dispersions at T, T, and K. 


On the Go Wo level of accuracy the CBM in bulk MoS2(T; point) 
is ~0.4 eV lower in energy than the CBE at K. and the VBM 
at [is favored by ~0.5 eV energy difference with respect to the 
band edge at K. In bulk MoS;, the valence band splitting at K 
amounts to roughly 240 meV. 

In contrast to bulk and multi-layer MoS», the main attribute 
of the single-layer MoS. band structure is the direct fun- 
damental band gap located at K since both, the VBM and 
the CBM, are at K. This direct band gap at K has been 
clearly demonstrated by several preceding works on single- 
layer MoS»(140] and explains the 


significant increase of photoluminescence yield in single-layer 


MoS»[I0]. The other important feature in the band-structure 
of IL-MoS, is the splitting of the valence band maximum 


Ky, — Ky2 which amounts to ~150 meV. Since inter-layer in- 
teractions are absent, this splitting must have a origin different 
from the splitting in bulk. Indeed, due to missing inversion 
symmetry, the spin-degeneracy of the bands is lifted, resulting 
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in a particularly large spin-orbit splitting at K[146][143]. This 
splitting explains the doublet structure of the strong peak in the 
absorption spectrum of 1L-MoS, [10]. 

The CBM in 1L-MoS; is also located at K but the splitting 
due to SOC is one order of magnitude weaker than the splitting 
of the VBM. Its absolute value is strongly affected by the ex- 
change correlation functional used in the calculations [147] as 
will be discussed later. Both, the valence and conduction bands 
exhibit nearly parabolic dispersion at this point, which explains 
the small effective charge carrier masses and indicates promis- 
ing conductivity. However, compared to bulk MoS, the va- 
lence band at T is considerably flattened. This flattening results 
in a much higher effective hole mass in 1L-MoS», which was 
also observed in Angle-Resolved Photoemission Spectroscopy 
(ARPES) experiments [148]. A second local conduction band 
minimum close in energy is observed at T. The relative energy 
positions of the states K, and T, and the location of the VBM 
(either K or I) determine, whether the material is a direct or 
an indirect semiconductor. We observed that in 1L-MoS» the 
Ke — T, energy difference is very sensitive to (i) the structural 
optimization, (ii) the applied in-plane strain, and (iii) the GW 
accuracy (around 0.05 eV). We discuss these issues in more de- 
tail later in this section. 

We now describe the changes stemming from quasiparticle 
corrections in the band-structures of bulk, single-, and double- 
layer MoS. The most notable change is the sizable increase of 
the band gap on the level of the GW method. Also the valence 
band width increases slightly. Note that for the single-layer, 
the VBM at T is shifted downwards as compared to the VBM 
at K. The conduction band is even more profoundly affected. 
The upshift of the CBM at K, is larger than the one of the sec- 
ondary CBM T.. In the single-layer, this results in the lowering 
of the T. energy relative to K, and thus in the reduction of the 


K.-T, energy difference by 60% compared to LDA turning the 
material almost indirect on the GoWo level. In bulk MoS» the 
CBM T. is lower in energy than Ke (due to inter-layer interac- 
tion) already on the DFT level. On the Go Wo level, this trend is 
amplified even more. In both, double-layer and bulk MoS;, the 
GW corrections 

Besides that, one has also to consider that the GW correction 
to the band structures slightly depends on the number of layers 
in multi-layer systems. We find that the band gap correction at 
K decreases in double-layer and bulk with respect to the case 
of single-layer MoS». The larger number of layers is associated 
with an increase of the dielectric screening, which results in a 
smaller correction |149}, [140]. 

In order to better understand, the origin of the differences be- 
tween single-, double-, and bulk MoS, band structures, we have 
summarized the orbital composition of the highest valence and 
lowest conduction bands at points of interest in the Brillouin 
zone in Table [5] In single-layer MoS», the S-p and Mo-d or- 
bitals dominate the composition of the wave functions, with a 
minor contribution of the s orbitals. The conduction band edge 
at K is mainly a Mo-d; (86 96) and the remaining part shared 
between the S-p,, and Mo-s orbital. The valence states at K,ı 
and K, are predominantly Mo-d,, (80 96) and 20 % of S- p, 
without any contribution from s orbitals. The wave function 
at the local minimum at T. has a more complex composition, 
typical for points of low symmetry, as summarized in Table 
These GW findings qualitatively reproduce previous DFT- 
PBE results (e.g., Fig. 4 in Ref. [28], Fig. 5 in Ref. [32] 
as well as the tight-binding (TB) model of Liu and co-workers 
(150][I51]. The latter, however, suggest also a significant con- 
tribution of the Mo-d,2_,2 states to the valence band edge at K, 
which we contribute to a deficiency of the TB model using only 
three Mo-d bands. The composition of the I' and K, states in 
bulk MoS» is very similar to the single-layer values. The bulk 
K,; and Kj» states however, are now predominantly composed 
of Mo-d,2_,2 and the valence band states at I', change the weight 
of the orbital p. of sulphur atoms. The latter increase is related 
to the bonding between sulphur atoms of different layers, which 
produces the interlayer coupling [152]. 

A correct description of the electronic properties requires (1) 
the inclusion of Molybdenum semi-cores states (4s and 4p or- 
bitals) in the basis set, (11) a plane wave cutoff of 350 eV, (iii) at 
least a 12x 12x3 (12x 12x 1) I-centered k mesh for bulk (1L) 
MoS), and (iv) the explicit inclusion of the spin-orbit interac- 
tion [153]. We interpolate the band structure to a finer grid us- 
ing the WANNIER90 code and the VASP2-WANNIER90 
interface [155]. With respect to GW calculations, it is impor- 
tant to mention that (1) solely including valence electrons leads 
to an erroneous wave-vector dependence of the GW correction 
[145], Gi) the convergence with respect to virtual states when 
calculating W is particularly slow for 1L-MoS»[156], and (iii) 
the default value for the number of quasiparticle energies that 
are calculated and updated in the scGW VASP calculation must 
be substantially increased. While taking into account « 500 
virtual states is sufficient to converge the quasiparticle gaps of 
bulk MoS» within 20 meV, more than 1000 bands are required 
for 1L-MoS, gaps to be stable within 40 meV. The logarithmic 
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Figure 11: Convergence of the direct and indirect quasiparticle energy gap of 
1L-MoS> calculated within the GoWo approximation with respect to the total 
number of bands (occupied + virtual) taken into account for the calculation of 
the screening (W). The solid lines represent logarithmic fits and serve as guides 
to the eye. 


scaling of the direct and indirect gap in 1L-MoS» with respect 
to the number of bands included in the calculation is shown in 
Fig. Concerning the number of quasiparticle energies that 
have been updated in the scGW calculations (NBANDSGW 
parameter in VASP), it is emphasized that more than 200 are 
required to converge the quasiparticle gaps. In particular the 
conduction band extremum at point T strongly depends on this 
parameter. 


4.2. Dependence on the crystal structure 


The analysis of the preceding paragraphs underlines the im- 
portance of accurately calculating the energy difference be- 
tween the conduction band minima at K, and Te. This is a 
challenge for the different theoretical approaches mentioned be- 
fore, because these quantities also sensitively depend on the de- 
tails of the crystal structure. In order to discuss this, we focus 
on single-layer and bulk MoS» but the conclusions can be ex- 
tended to multi-layer MoS;. 

One source of controversy between the results of several cal- 
culations performed for 1L-MoS; could 
be the underlying crystal structure. In particular, the lattice 
constant, interlayer distance (relevant for multi-layer and bulk 
MOS»), and atomic positions defining the interatomic distance 
(Mo-S bond length and S-Mo-S bond angle) may significantly 
affect the energy gaps and band dispersion. The dependence 
of the MoS band structure on the details of the crystal struc- 
ture has not been addressed so far and will be elucidated in the 
following. 

Most calculations reported so far, have used the experimen- 
tal room temperature lattice constant of bulk MoS,, i. e., 
a = 3.16 Å, and 10-15 A vacuum along the c axis for the single- 
layer (1L) MoS, slab structure. However, less information is 
given about the choice of the origin of the unit cell (atomic po- 
sitions) and the z parameter. Unfortunately, according to Bron- 
sema et al. an inconsistency exists in literature regard- 
ing the choice of the atomic positions and the corresponding 


Table 5: Orbital composition of the wave functions at the points K, T and T, in the case of single-layer and bulk MoS2. 


Single-layer 


Atom Sulphur Molybdenum 

Orbital s Pry Pz s de dy do dy dy 
re - 0.54 - - - 0.46 - - - 
T, - - 0.23 | 0.02 - - 0.75 - - 
Ke - 0.09 - 0.05 - - 0.86 - - 
Ky - 0.20 - - - - - - 0.80 
Ky - 0.20 - - - - - - 0.80 
Te 0.03 0.22 0.06 - 0.54 - 0.12 - 0.01 
Bulk 

Atom Sulphur Molybdenum 

Orbital s Pxy Pz s dey dy dz dwz dwy 
I, - 0.53 - - - 0.47 - - - 
D, 0.07 - 0.30 | 0.03 - - 0.60  - - 
Ke - 0.09 0.05 - - - 0.86 - - 
Ky - 0.21 - - 0.79 - - - 

Ky - 0.18 - - 0.82 - - - 

T, 0.07 0.18 0.06 - 0.52 - 0.14 - 0.03 


z parameter. In fact, the latter determines the interatomic dis- 
tances and thus the S-Mo-S layer thickness. For this reason, 
the band structure of bulk and 1L-MoS; has been calculated by 
LDA+SOC and Go9Wo-4 SOC for some of the crystal structures 
summarized in Tab. |1| The GoWo--SOC results are depicted in 
Fig. 

For bulk MoS;, small variations of the z parameter (com- 
pare red and blue lines in Fig. and the lattice constants a 
and c (compare the red and green lines in Fig. have only 
little influence on the band-structure, since there is anyway a 
strong inter-layer interaction that leads to a splitting of the va- 
lence and conduction bands. The situation is quite different for 
the single-layer: as illustrated in Fig. BO The significant role 
of the internal z parameter that is defining the interatomic dis- 
tances (Mo-S bond length and S-Mo-S bond angle) is revealed. 
With increasing z from 0.621 to 0.629, the Mo-S bond length 
is reduced from 2.42 Å to 2.35 Å, respectively. This favors 
hybridization between the Mo-d and S-p states that comprise 
the highest occupied and lowest unoccupied bands and there- 
fore the band dispersion (band width) increases. As a conse- 
quence, the CBE at T and the VBE at T are pushed to higher 
energies and the CBE at K becomes the CBM giving rise to a 
direct fundamental energy gap. It is strongly emphasized that 
an improper choice of atomic positions and corresponding z pa- 
rameter can fortuitously yield a direct fundamental band gap 
in IL-MoS, and may partly explain the inconsistency among 


Go Wo band structures|140 reported so far. 


Another source of discrepancies between single-layer 
MoS, calculations might be related to the relaxation of the 
atomic positions. In Fig. Ba), the GoWo band structures of 
1L-MoS, calculated using the experimental bulk unit cell lat- 
tice constant (a=3.16 Å) without and with LDA-relaxed atomic 
positions in the single-layer are compared. Due to the overbind- 
ing in DFT-LDA, the Mo-S bond length gets reduced with re- 
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Figure 12: (Color online) Band structures of bulk (a) and single-layer 
MoS, (b)calculated within the Go Wo approximation explicitly including SOC 
for different experimental crystal structure parameters. The valence band ex- 
tremum at K is aligned at zero energy. The abbreviations ACB39 (3.161 À, 
c = 12.295 A, z = 0.6275), JCG15 (3.14 A, c = 12.327 A, z = 0.621), ZAC540 
(3.16 A, c = 12.29 A, z = 0.629) refer to the crystal structure details published 
in Ref. [71], Ref. (72), and Ref. [157], respectively. 


laxation, which again strengthens the Mo-d-S-p hybridization 
resulting in an increase of the band dispersion along PK. This 
results in a raise of the VBM at I accompanied by an increase 
of the conduction band T valley energy and the stabilization of 
the CBM at K. Consequently, GgWo+SOC yields the correct 
direct gap band structure for 1L-MoSg, if the atomic forces are 
minimized on the LDA level. As can be seen from Fig. Bta, 
the direct gap at K is reduced in the position relaxed case by 
~90 meV. While the CBM energy at T is not affected, further 
reduction of the direct gap at K by ~40 meV is obtained by us- 
ing the lattice constant of the optB86b-VdW fully relaxed bulk 
unit cell (a=3.164 Å) and LDA-relaxed atomic positions in the 
single-layer [not shown in Fig. [I3[)]. 

A final test on the level of Go Wo for the influence of the struc- 
tural details on the band structure of 1L-MoS? was performed 
with a fully optB86b-VdW optimized single-layer structure, 
i. e., the in-plane lattice constant a = 3.162 À and LDA-relaxed 
atomic positions. Since the VdW interactions are not relevant 
in the single-layer, the obtained structure is very close to the 
experimental bulk one. Thus the Go Wo+SOC band structure re- 
sembles that one calculated without any atomic position relax- 
ation [indirect gap, not shown in Fig. Bta]. From this analysis 
we conclude, that the location of the valence and conduction 
band extrema at I and K are very sensitive to the relaxation of 
the atomic positions and if the atomic positions in the single- 
layer are relaxed within DFT-LDA, the CBM at K is stabilized 
with respect the CBE at T. 

The finding of an direct gap 1L-MoS, on the GoW level 
seemed to be controversial to the results obtained by Shi 
et al. (144], who performed an analogous comparison of GoWo 
and scGW calculations for 1L-MoS» as presented in Fig. mfb) 
and concluded that single-shot Go Wo is insufficient in describ- 
ing 1L-MoS,. As shown here, omitting the relaxation of the 
atoms in the single-layer structure constructed from the exper- 
imental bulk structure leads indeed to an incorrect description 
of IL-MoS» on the GoWo level (indirect band gap). The scGW 
calculation cures this problem, but further increases the direct 
gap at K. The tendency of scGW to overestimate semiconduc- 
tor band gaps is known and thus one must assume that the 
scGW direct gap of 1L-MoSz is too large. 


4.3. Performance of different methodologies 

After the discussion of the relation between structural and 
electronic properties, we focus on how the results depend on the 
XC functionals. The results of this analysis for bulk, double-, 
and single-layer MoS, are summarized in Table [6]and Fig. 
Note, the fully optimized structures using the van der Waals 
functional as previously described were used in these calcula- 
tions. 

On all levels of theory the band structure of bulk 
MoS shown in Fig. [14[a) and (b) corresponds to an indirect 
semiconductor with the VBM located at I' and the CBM at T. 

When comparing the numbers given in Table le] LDA and 
PBEsol underestimate the T, — T, transition (indirect gap of ~ 
0.85 eV) are found to be in agreement with previous calcula- 
tions [142] [141]. Compared to LDA, the inclusion of local ex- 
change as provided by the MBJ potential mainly affects the T, 
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Figure 13: (Color online) Band structure of IL-MoS» calculated with the 
experimental bulk lattice constant of a=3.1602 Å. In (a) the Go Wo band struc- 
ture obtained by omitting the atomic force minimization in the single-layer is 
compared to the corresponding results with atomic position relaxation on LDA 
level. In (b) the Gp Wo+SOC and the scGW-- SOC approach are compared. The 
VBM at K is set at zero energy. The Fermi level is indicated by the dashed 
horizontal line. 


and I, energies resulting in a larger indirect T. — T, gap of 1.15 
eV. However, the band dispersion along K towards T is reduced 
resulting in almost energetically balanced CBEs at K and T, 
i. e., |K. — T.| is only 60 meV. Therefore the difference between 
indirect Te — T, and the direct Ke — K, gap is less pronounced 
than in LDA. Further improvement is achieved by taking into 
account non-local exchange using the HSEsol functional that 
increases the indirect as well as the direct Ke — K,, gap. The 
|K.—T,| difference in HSEsol is roughly 140 meV, significantly 
less compared to LDA (~ 240 meV). Taking into account the 
dielectric screening in the GoWo approach strongly enhances 
the |K; — T.| difference to ~ 390 meV. This results in an indirect 
for bulk MoS, of 1.24 eV, in good agreement with experiment 
(1.2-1.3 eV) as well as previous calculations. The 
outstanding agreement of both, the indirect (1.24 eV) and direct 
(2.08eV) bulk gaps, with the values obtained by an all-electron 
GW code verify the accuracy of the present PAW results. 
A higher level of accuracy is reached by performing scGW 
calculations. Compared to the Go Wo band structure, the |K.— T.| 
difference is again significantly reduced (to ~ 160 meV), since 
T, is pushed up in energy almost back to the HSEsol position. 
The fundamental indirect scGW gap is 1.39 eV and slightly 
overestimated compared to experiment, which is due neglect- 
ing the attractive electron-hole interactions via Vertex correc- 
tion in GW. In the K-T region, the scGW band structure 
resembles the HSEsol, whereas we observe remarkable differ- 
ences in the I-M-K range. The GoWo and scGW results are 
consistent to previous calculations given in Tab. [6] within 
the uncertainties originating from computational aspects. 
Single-layer MoS% is described as a semiconductor with a di- 
rect gap at K on all levels of theory beyond standard DFT-LDA, 


Table 6: Direct band gaps and interband transitions in MoS (in eV) as well as the energy difference between the two lowest conduction band extrema Ke — Te 
calculated on different levels of theory explicitly including SOC in comparison to available literature data listed in brackets. 


(eV) Ke m Ky Ke = Ky Ke = D, Te 7 Ky Te E D, Ke zi Te Te m LL, 
Bulk 

LDA 1.64 1.86 1.07 1.40 0.83 -0.24 2.08 
LDA (1.80) (0.81) 

PBEsol 1.65 1.87 1.10 1.42 0.87 -0.23 2.11 
PBE (0.87) 

PBE (1.58) (0.86) 

MBJ 1.62 1.82 1.21 1.56 1.15 -0.06 2.27 
HSEsol 2.10 2.36 1.58 1.96 1.44 -0.14 2.84 
HSE (2.16) (1.48) 

GoWo 2.08 2.32 1.63 1.69 1.24 -0.39 2.53 
GoWo (2.07) (1.23) 

GoWo (2.00) (1.30) 

scGW 2.17 2.41 1.59 2.02 1.44 -0.16 2.88 
scGW (2.099) (2.337) (1.287) 

EXPT. (1.78) (1.29) 

EXPT. (1.95) (1.20) 

Single layer 

LDA 1.62 1.77 1.61 1.92 1.91 0.30 2.74 
PBEsol 1.65 1.79 1.69 1.89 1.93 0.24 2.78 
PBE (1.75) 

PBE (1.60) 

HSEsol 2.09 2.28 2.23 2.45 2.59 0.36 3.63 
HSE 2.06 2.25 2.13 2.51 2.58 0.45 3.60 
HSE (2.05) 

GoWo 2.45 2.60 2.61 2.59 2.74 0.14 3.60 
GoWoll41) (2.82) 

GoWo (2.97) (3.26) 

scGW 2.72 2.87 2.90 2.98 3.16 0.26 4.29 
scGWo (w/o SOC) (2.78) 

scGW (2.759) (2.905) 

EXPT. (1.88) (2.05) (1.6) 
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provided that the crystal structure is fully relaxed as stressed in 
Sec. Standard DFT (LDA and PBEsol) severely underesti- 
mates the direct gap of the single-layer structure. Besides that, 
LDA wrongly sets the VBM at I at slightly higher energy than 
the VBE at K. It is important to emphasize that the underes- 
timated K,ı — I, energy difference is observed if the single- 
layer slab is constructed from the optB86b-VdW relaxed bulk 
structure (in-plane lattice constant a=3.164 A), but not in case 
of the experimental bulk lattice constant a=3.16 A. Therefrom 
we conclude that the relative positions of the K,; and T, en- 
ergies are strongly dependent on the in-plane lattice constant. 
Hence is imperative the investigation of strain effects on the 
1L-MoS, band structure presented later in this section. 

Employing the HSEsol functional to 1L-MoS» shifts the 
conduction bands almost uniformly upwards in energy com- 
pared to DFI-LDA resulting in a rather constant |K; — Tel, 
i. e., ~360 meV and ~300 meV, in HSEsol and LDA, respec- 
tively. Band dispersions in the valence band region are in- 
creased within the HSEsol description and the VBM splitting 
at K is enhanced to ~ 200 meV compared to the LDA value 
of ~ 150 meV. The VBM at K,, is stabilized by roughly 200 
meV compared to I, in HSEsol calculations. Concerning GW 
approaches, the subtle changes of the band dispersions between 
T and K result in a significant change of the K, — T, energy dif- 
ference: Go Wo stabilizes the CBM at K by 130 meV, whereas 
scGW enhances this energy difference by a factor of two. Note 
that the energy differences strongly depend on the total num- 
ber of bands (NBANDS) included in the GW calculations (see 
Fig. and the convergence with NBANDS itself is influ- 
enced by the amount of vacuum included in the single layer 
MoS) cell (20 A in the present case). This means that using a 
larger amount of vacuum requires an increase of the NBANDS 
parameter as well. For this reason, the comparison between the 
present results and previously reported values, as summarized 
in Tab. [6] is difficult. The values listed in Tab. [6]refer to GW 
calculations with NBANDS=512. Increasing NBANDS from 
512 to 1920 reduces the direct gap at K by roughly 80 meV, 
the indirect gap by 60 meV, but the Ke — T. energy difference 
increases by 40 meV. 

Analogous to bulk MoS», including non-local exchange by 
HSEsol increases the Ke — K,; gap (2.09 eV) considerably. 
The calculated Ke — K,; quasiparticle GoWo gap amounts to 
2.45 eV, which is smaller by 0.3-0.5 eV to reported values. 
This difference is attributed to structural and 
computational details: A Go Wo» calculation performed with the 
experimental crystal structure and a reduced k mesh of 8x8x2 
yields 2.86 eV. Liang et al.reported a direct band gap of 1L- 
MoS, of 2.75 eV, [160] which was obtained by GoWo calcu- 
lations taking into account a Coulomb interaction truncation 
to avoid spurious interlayer interaction between the periodi- 
cally repeated monolayers, but using the generalized plasmon- 
pole model (GPP) for the dynamical screening and omitting 
SOC. The issues of the Coulomb interaction truncation, k-point 
sampling, and vacuum layer thickness were also addressed by 
Hüser et al.[161], who argued that the band gap values con- 
verged with respect to k-point sampling and slab distance are 
rougly 0.4 eV too small compared to the free standing mono- 
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layer (including Coulomb truncation). Once again this re- 
flects the difficulty to achieve accurate results and explains the 
plethora of band gap data in literature. 

Compared to the single-shot Go Wo result for the direct Ke — 
K,ı gap, the scGW further increases the gap to 2.72 eV in 
agreement with previously reported values as listed in Table [6] 
The slow convergence of the IL-MoS; GW band gaps with 
the NBANDS parameter was only recently stressed and 
might explain the larger values previously reported. 

At this point it is important to recall that quasiparticle gaps 
are single-particle gaps. Their overestimation by roughly 0.5 
eV compared to some experiment is explained by the missing 
electron-hole interactions (excitonic effects), which are strong 
in 2D materials due to confinement and lead to the formation 
of bound electron-hole pairs. These bound excitons reduce the 
direct band gap by their binding energy and define the optical 
gap, which is experimentally accessed by optical measurements 
such as photoluminescence or photoconductivity. Excitonic ef- 
fects are addressed in the next section. 

To conclude, it is evident based on the above discussion that 
the different levels of accuracy and/or complexity applied in the 
GW methods substantially alter the results. The non-local ex- 
change and dynamical screening are inevitable for an accurate 
description of the electronic properties of MoS». Based on our 
calculations of single-layer MoS, a reasonable estimate for the 
direct band gap is 2.4+0.2eV and the spin-orbit splitting of the 
valence band edge at K is 150-160 meV. The energy difference 
between the two valence band extrema in 1L-MoS, is much 
smaller than in bulk MoS, and very sensitive to the in-plane 
lattice constant. As put forward by Kuc and Heine[147], the 
estimations of the weak spin-orbit splitting of the conduction 
band edge is strongly dependent on the XC functional used in 
the DFT calculation and a better description by methods beyond 
ground state DFT is required. From our GW calculations, we 
deduce a conduction band splitting of 10 meV, which however 
falls within the estimated uncertainty range. 

A final remark concerns the starting point of the scGW cal- 
culations. One should keep in mind that the scGW result can be 
influenced by the wave functions (orbitals) used for calculating 
G and W as pointed out in Ref. [162]. Thus, in the complex- 
ity of the GW method one can go a step forward by applying a 
full self-consistent procedure in the Green's function G and the 
screened interaction W, consisting in the iterative solution of 
the Dyson equation. However, the extraordinarily cumbersome 
calculations required for this procedure restrict the application 
of this approach to small systems, such as binary molecules as 
N» or CO [162]. Up to now this scheme has not been applied 
to single-layer MoS» and we think that its implementation for 
layered materials is still far. 


4.4. Strain effects in single-layer MoS; 


The ideal scenario of free-standing 2D layers as considered 
in most theoretical simulations is hardly fulfilled in reality. In 
the course of experiments or device fabrication with 2D ma- 
terials, it is important to consider strain resulting e.g. from 
the lattice constant mismatch between the substrate and the 


2D layer. Equally important in this context is the interac- 
tion of the 2D material with the substrate as shown in Ref. 
[148]. Therein, the ARPES scans of exfoliated single-layer 
MoS, compared to those of chemical vapor deposition grown 
single-layer MoS, on silicon revealed that the presence of sub- 
strate alone is sufficient to modify the MoS, band structure. In 
particular, the MoS>-substrate interactions are responsible for 
the pronounced flattening of the VBM at I of MoS2 on silicon. 

In addition, recent experiments have demonstrated that appli- 
cation of tensile strain changes the gap from direct to indirect 
(30]. In particular, the MoS, flake deposited on a flexible sub- 
strate which is subsequently deformed in a controllable manner, 
experiences uniaxial tensile strain up to 2.2 %. The photolumi- 
nescence spectra of these samples show a clear transformation 
of the band character, and an associated reduction of the inte- 
grated intensity of the optical signal. 

The sensitivity of TMDs band structure on the lattice con- 
stant opens the possibility to modify the band gap and thus the 
optical properties in a controlled way by external strain. This is- 
sue has been theoretically addressed either through LDA/GGA 
calculations [31], or the GW method [144]. The 
effect of hydrostatic pressure on the vibrational, electronic, and 
optical properties of bulk, multi-, and single layer MoS. was 
investigated by Nayak et al. by combining various ex- 
periments (high resolution transmission electron microscopy, 
electrical resistance measurements, laser Raman spectroscopy, 
synchrotron X-ray diffraction experiments under high-pressure) 
with DFT calculations. Interestingly, while the direct bulk band 
gap decreases with increasing pressure, the direct band gap of 
1L-MoS, increases by 11.7% up to ~ 12 GPa before it is re- 
duced. Thus the pressure induced electronic transition from the 
semiconducting to a semimetallic state occurs at much larger 
pressures in the latter. [34] 

Being aware of the importance of substrate interactions, we 
investigated the strain effects on the electronic properties of 1L- 
MoS, within the GW approach and the model of free-standing 
2D layers. Biaxial tensile strain has been realized by increasing 
the in-plane lattice constant of 1L-MoS». The band structures 
of the strained materials were calculated with relaxed atomic 
positions and are shown in Fig. The direct Ke — K,, gap and 
interband transitions as a function of strain deduced from these 
band structures are collected in Tab. 

With increasing in-plane lattice constant (biaxial tensile 
strain) the bond distances within the xy plane of the Mo-S- 
Mo sheets are changed, but also in the perpendicular direction 
through the relaxation along z. According to Tab. |5| the valence 
states at K (VBM) are mainly composed of S-p,, and Mo-d,, 
orbitals, whereas the valence states at [ have predominantly Mo 
d.» character. Concerning the CBM at K, the states are mainly 
Mo d; orbitals and the conduction band states around T have 
predominantly Mo d» ,» character. By changing the S-Mo-S 
bond lengths and angles due to tensile strain, the overlap of the 
Mo dz with the S p,, is reduced, whereas the coupling between 
the Mo d,, and S p,, is increased.[32] As a consequence, the 
TI, energy raises with respect to K, and the K, energy decreases 
compared to Te. Concomitantly, the Ke — K,, gap decreases, 
but not as fast as the K. — I, gap, which results in a transition 
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Figure 16: (Color online) The K; — K,; gap and the Ke — T, gap as a function 
of uniaxial [100] (dashed lines) and biaxial (full lines) tensile strain calculated 
with GoWo--SOC. The obtained transition strain is marked by the vertical full 
(dashed) lines in green for biaxial and uniaxial strain, respectively. The un- 
strained material refers to the experimental lattice constant a=3.1602 A. The 
experimental points from Ref. are represented by black dots. Note, the 
experimental data have been constantly shifted upwards in energy, since the 
electron-hole interaction is not included in the calculations. 


to an indirect 1L-MoS; as illustrated in Fig. Also evi- 
dent in Fig. [16] is the linear dependence of the band gaps on 
biaxial tensile strain. GoW +SOC suggests a transition strain 
of ~1.6% (i. e., 3.21 A), which is in good agreement with the 
value obtained by recent GW) calculations (~1.5%) [144]. Thus 
iterating the QP energies and one-electron wave functions in 
G only, seems to change the linear decrease of the band gaps 
only marginally. The scGW approach though, increases the di- 
rect Ke — K,, gap rather constantly (by ~ 0.2 eV) and signifi- 
cantly affects the band dispersions. The former results in a rigid 
shift of the direct gap dependence with strain, whereas the latter 
changes the K; — T, and K,, — T, energy differences (compare 
data summarized in Tab. p. Consequently, a much lower strain 
of ~0.7% (i. e., 3.18 A) for the direct-to-indirect transition in 
1L-MoS), is obtained. 

At this point, the present GW results are also compared 
to conclusions drawn from previous DFT-PBE calculations. 
Scalise et al.[27] proposed a stain-induced semiconductor to 
metal transition in 1L-MoS, on the basis of DFT-PBE calcu- 
lations omitting SOC. When the slope of the linear fit to their 
band gap data as a function of biaxial tensile strain (Fig. 2(b) 
in Ref. (27]) is compared to that deduced from the present GW 
results, the obtained difference is roughly 30%., i. e., the DFT- 
PBE slope is 30% smaller. However, in both cases the closure 
of the indirect band gap with increasing tensile strain is esti- 
mated at (10+1)% tensile strain suggesting that the trend of the 
electronic properties as a function of strain is well reproduced 
by standard-DFT. 

To compare with recent photoluminescence data [30], we 
also calculated the band structure of 1L-MoS» as a function of 
uniaxial tensile strain along the [100] direction with the Go Wo 
approach.These results are summarized in Tab. Since the 
electron-hole interaction is not included in the calculations at 


Figure 15: (Color online) The band structure of 1L-MoS» as a function of biaxial (left) and uniaxial (right) tensile strain calculated with Go Wo+SOC is depicted. 
The VBM at K is aligned at zero energy. The unstrained material refers to the experimental lattice constant a=3.1602 A. 


Table 7: The Ke — K, and Ke — T, gap as well as the Ke — Te energy difference in units of (eV) as a function of biaxial tensile strain. 


strain (%) GoWo+SOC scGW+SOC 

Ke S Ky Ke m D, (eV) Ke m T. Ke a Ky (eV) Ke = D, (eV) Ke = Te 
0 2.51 2.71 0.08 2.76 2.85 0.39 
1 2.37 2.44 0.24 2.61 2.58 0.55 
2 2.23 2.18 0.41 2.47 2.31 0.72 
3 2.11 1.93 0.55 2.34 2.06 0.86 
4 1.98 1.69 0.69 2.21 1.82 0.99 
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this stage, the experimental data have been rigidly shifted in en- 
ergy for better comparison. Thereby it is assumed that the band 
gap renormalization due to excitonic effects does not change for 
small strains as applied in this case. 

The effect of the uniaxial strain on the band structure as 
shown in Fig. [I5]is similar to that of biaxial strain. Again the 
K, energy decreases, whereas the I, energy increases with in- 
creasing uniaxial tensile strain giving rise to the direct-indirect 
gap transition, when I’, becomes the VBM. Compared to biax- 
ial strain, the slope of the linear dependence of the direct and 
indirect energy gaps is significantly smaller in the uniaxial case, 
i. e., roughly by a factor of two. As a consequence, the transi- 
tion between direct and indirect single-layer MoS» occurs at a 
larger strain of 3.3% (equivalent to a=3.26A), which is roughly 
twice as high as in the biaxial case. Since the strain induced 
transition to an indirect gap significantly reduces the photo- 
luminescence yield, understanding of the strain effects on the 
opto-electronic properties of single-layer MoS» is particularly 
relevant for technological applications. 


4.5. Tight-binding modelling of single-layer MoS 

Tight-binding calculations can give further insight on the 
electronic properties of single-layer MoS». The tight-binding 
method expands the wave functions in terms of an atomic or- 
bital basis, thus giving a simple and intuitive physical picture 
of the electronic bands. The atomic orbital weight of each band 
state is directly accessible and changes in the band structure can 
be attributed to the change of a single tight-binding parameter. 

Fig. shows the change of the conduction and valence 
bands with variation of the atomic positions. However, it does 
not tell which orbitals are responsible for such variations. The 
parametrization given in Ref. was used, but for clarity, 
we only consider d-orbitals for Mo and p-orbitals for S. In- 
teratomic interactions up to the second nearest neighbors were 
taken into account. Figure[17]illustrates the band structure alter- 
ing the hopping parameter V while keeping fixed the remaining 
parameters. The subindices denote the kind of orbitals (p or d) 
and the symmetry of the bond, ct, and delta (see details in 
(1641[165]). The size of the red circles indicates the weight of 
the do orbital and that of the blue circles the weight of the dz 
orbitals. We have increased each parameter +10%. 

Focusing on nearest neighbor interaction between Mo-S, the 
parameters V4, and Vj; carry the effect of a vertical displace- 
ment of the S atoms. In the ab initio calculations shown in Fig. 
[I2]we have seen the total effect of the displacement on the band 
structure, but with the tight-binding formalism we can assign 
interactions to bands. The change in V,q, modifies the rela- 
tive positions of the two local minima of the conduction band, 
K, and T,, producing almost a rigid shift of the valence band. 
In the case of the hopping parameter V,4;, the effect is rather 
different, giving a rigid shift of the conduction band at K, and 
Te. For the valence band, the dispersion at I', is strongly shifted 
to higher energies for +10% and to lower energies for —1096 
Vpan. The combination of these two weights is equivalent to 
displacing vertically the S atoms. In the case of the hopping 
parameters related with the d orbitals of Mo atoms, the results 
are more complex. This can be seen for the cases of Vaas and 


27 


Figure 17: 
parameters (indicated in each graph). For each graph, we vary the indicated 
parameter +10 96. The red circles are associated with the do orbitals and the 
blue circles with the d? orbitals. Their size is related to the orbital weight. 


Band structures for different values of the tight-binding hopping 


Vaas in Fig. Only for the Via; case a minor influence on the 
band structure is obtained. The latter three parameters are also 
important in the change of the lattice parameter. The analysis in 
terms of the tight-binding model establishes the rather complex 
nature of covalent bonding formed by p — d orbitals in TMDs. 
Note that a proper description of the conduction band would 
require a tight-binding model with more parameters to capture 
the interaction with farer neighbors and to enlarge the basis of 
atomic orbitals [166]. 


4.6. Effective charge carrier masses 


Electron and hole carrier mobility is inversely proportional to 
their effective masses. The strain effects on effective mass are 
an important issue from the technological point of view. The ef- 
fective electron and hole masses of bulk MoS, as a function of 
strain were studied on the DFT-HSE level by Peelaers and Van 
de Walle and on DFT-LDA level by Scalise et al..[28] 
For the sake of completeness, the values for the bulk effective 
masses given in Ref. are briefly summarized here. As 
pointed out in Ref. and evident from Tab. III of Ref. 
(123], the valleys as well as edges at K exhibit rather isotropic 
effective masses concerning both longitudinal and transversal 
directions. This justifies averaging over both directions as sug- 
gested in Ref. [141]. However, the longitudinal and transver- 
sal masses are quite different at the T point (the same as A 
in Ref. and Emin in Ref. ([143]) and the I point, which 
are the CBM and VBM in bulk MoS;, respectively. The value 
given in Ref. for the transversal hole mass (i. e., I [X] 
in Ref. (123]) is 0.62mo, where mo is the electron rest mass. 
This number is in good agreement with the present result of 
0.64mo, evaluated by a parabolic fit of the Wannier interpolated 
GoWo+SOC band structure, calculated with the experimental 
structure data according to E = IE. in the I'M direction (equal 
to TK). For the optB86b-VdW relaxed bulk, we obtain 0.69mo 
and 0.70mp on the GgWo+SOC and scGW+SOC level, respec- 
tively, indicating that the effective masses are well described on 
the GoWo--SOC level. The longitudinal hole mass given in Ref. 


Table 8: The Ke — K,; and Ke — T, gap as well as the Ke — T; energy difference in units of (eV) as a function of uniaxial [100] tensile strain. 


strain (76) GoWo4 SOC scGW+SOC 

Ke = Ky Ke = D, (eV) Ke m Te Ke = Ky (eV) Ke = T, (eV) Ke = Te 
0 2.51 2.71 0.08 2.85 2.76 0.39 
1 2.43 2.56 0.19 2.58 2.56 0.53 
2 2.37 2.44 0.27 2.47 2.49 0.60 
3 2.29 2.31 0.36 2.33 2.42 0.69 
4 2.21 2.17 0.46 2.20 2.35 0.77 


| 123] is 0.80mo reflecting the anisotropy between transveral and 
longitudinal hole masses. The corresponding estimate from the 
GoWo+SOC band structure of the optB86b-VdW relaxed bulk 
is 1.05mo and 1.037 for the experimental bulk structure. This 
20% overestimation can be explained by the neglect of the spin- 
orbit interaction in the DFT-HSE calculations of Ref. [123]. 


The value given in Ref. for the transversal electron 
mass (i. e, Amin[A] in Ref. (123]) is 0.53mọ that is close to 
the present calculations, which yield 0.58mo at T (through fit- 
ting of E(k) along TK or TI). When evaluating the effective 
masses for both, electron and hole, at the K point, their longi- 
tudinal and transversal component are strongly isotropic. On 
the GoWo--SOC level, an average hole mass m, of 0.40mo and 
an average electron mass Me of 0.63mo is obtained. Perform- 
ing GoWo4 SOC (scGW+SOC) calculations with the optB86b- 
VdW optimized bulk structure yields m, = 0.39m (mj = 
0.40mo) and m, = 0.52mo (me = 0.52mo), for the average hole 
and electron masses, respectively. For comparison, the corre- 
sponding values reported in Ref. are 0.45mo and 0.46mo 
for m; and me, respectively. Despite the different theoretical ap- 
proaches, i. e., DFT with the HSE functional omitting spin-orbit 
coupling in Ref. compared to GW calculations including 
SOC in the present work, the overall agreement between the 
results is good. 


In the following, we focus on single-layer MoS» and present 
the effective electron and hole masses at the K point for the 
longitudinal and transversal directions in Tab. p] For compar- 
ison with available literature data, the average of the longitu- 
dinal and transversal component of the electron (me) as well 
as hole mass (715) are included too. Note that both, Go Wo and 
scGW calculations were performed with the optB86b-VdW op- 
timized in-plane lattice constant and LDA-relaxed atomic posi- 
tions. The values corresponding to the single-layer constructed 
from the experimental bulk structure are given for comparison 
in brackets in Tab. D] As emphasized by Shi et al. [144], the 
effective masses are sensitive to strain (i. e., the in-plane lat- 
tice constant), spin-orbit interaction, the GW accuracy as well 
as the convergence criteria of the GW calculations, in particu- 
lar the k point sampling. It is difficult to disentangle these de- 
pendencies and the values summarized in Tab. D]exhibit some 
spread. However they are consistent within the order of magni- 
tude. The average values obtained in the present work for the 
effective electron mass me range from 0.35-0.40mop, whereas av- 
erage effective hole masses m, between 0.43 and 0.49mọ have 
been estimated. 
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In contrast to the electron masses, the anisotropy between 
longitudinal and transversal hole masses is strongly pro- 
nounced, thus yielding a larger uncertainty range for the av- 
erage value. This may be explained by the significance of spin- 
orbit interactions: SOC affects the valence band dispersion at K 
more strongly than the CBM [144]. A common observation is, 
that the electron mass is slightly smaller than the effective hole 
mass. As shown in the previous section, tensile strain shifts the 
valence band maximum to the I point. At this point the band 
dispersion is significantly smaller than at the K point resulting 
in much larger effective hole masses (2-3mo estimated at 2% 
biaxial tensile strain) and therefore decreased hole mobilities. 


Table 9: Values of 1L-MoS> effective charge carrier masses in units of the electron rest mass mo evaluated from parabolic fits of the valence and conduction band 
edges. The subscripts refer to the longitudinal (/) and transversal (f) directions that are further specified by the points given in parenthesis. The first one denotes 
the location of the band extremum, whereas the second defines the direction from that point. As suggested in Ref.[141], the average effective masses (me and mp) 
determined from the longitudinal and transversal directions are also included. The values obtained from band structures calculated with the experimental in-plane 
lattice constant are given in brackets. 


Single-layer electron hole 
m(K-T) m,(K-M) m m(K-T) m(K-M) m, 

GoWo* SOC 0.33 0.36 0.35 0.40 0.50 0.45 

[0.36] [0.36] [0.36] [0.35] [0.50] [0.43] 
Go Wo+SOC [0.60] [0.54] 
HSE [0.37] [0.38] [0.38] [0.44] [0.48] [0.46] 
G,Wo+SOC 0.37 0.21 
scGW+SOC 0.39 0.42 0.40 0.42 0.56 0.49 

[0.33] [0.36] [0.35] [0.38] [0.50] [0.44] 
scGW+SOC [0.34] [0.35] [0.35] [0.46] [0.43] [0.44] 
scGWo 0.36 0.39 
scGWo [144 [0.32] [0.37] 
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5. Optical properties 


In semiconductors like MoS, electron-hole pairs interact via 
Coulomb attraction, forming excitons. Excitonic effects deter- 
mine the optical properties of MoS [167] [168] [169] [11|[10]. For 
example, experiments like photoluminescence and 
second harmonic generation are strongly influenced by 
excitonic effects. The most common excitonic effects are a red- 
shift in the optical gap (with respect to the quasiparticle gap) 
and, in some cases, a radical change in the optical spectra shape 
with respect to the independent particle spectra. This is in par- 
ticularly the case when bound excitons (absorption peaks below 
the onset of the continuum) are formed. It has been shown be- 
fore for hexagonal boron nitride (a prototypical wide-band gap 
layered material), that the anisotropic dielectric constant and 
the layered, quasi 2D confinement of excitons, leads to very 
strongly bound excitons[171] [149] [172]. Also in MoS;, there is 
a series of strongly bound excitons [145] [156] (albeit 
with comparatively lower binding energies). Their binding en- 
ergy depends on the number of layers (and the inter-layer dis- 
tance in the case of a periodic supercell calculation for single- 
layers). These excitons determine the shape of the optical ab- 
sorption spectra as will be explained in this section. 

We analyze the optical properties of MoS» multi-layers in 
the framework of the Bethe-Salpeter equation. In addition, we 
compare MoS, single-layer optical properties with other transi- 
tion metal dichalcogenides. Finally, we discuss results obtained 
applying empirical model Hamiltonians [174]. 


5.1. Bethe-Salpeter equation formalism 


The Bethe-Salpeter equation (BSE) formalism gives an accu- 
rate description of the electron-hole interaction [175]. BSE is 
based on many-body perturbation theory [176]. Starting from 
the eigenvalues and wave functions of the system, obtained by 
ab initio methods, BSE gives the dielectric function and the ex- 
citonic binding energy without introducing any additional pa- 


rameter |177]! 179]. BSE can be written as: 


(Ex — EWAN + 2, ErerlKenlevew) Avo: = OFAN. (15) 
k'v ^c mt 


The electronic excitations are expressed in the basis of 
electron-hole pairs, £,4. We assume vertical excitations at k, 
from a valence-band state with quasiparticle energy &y, to a 
conduction-band state with energy &x. AT. denote the expan- 
sion coefficients of the excitons and O* is the exciton energy. 
The interaction kernel Ken describes the screened Coulomb and 
the exchange interaction between electrons and holes, which 
includes local field effects [179]. In absence of electron- 
hole interaction Ken = 0. In this review, we consider only in- 
terband transitions. This is consistent with the experimental 
data, where the energy of optical excitations is always above 
the band gap value. Another important physical aspect is the 
omission of phonon-mediated transitions. They are important 
in indirect semiconductors, especially in the study of photolu- 
minescence [180]. We focus mainly on the optical absorption 
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spectra, where the weight of direct transitions is much higher 
than indirect transitions. 

The exciton wave function, expressed in the basis of the 
electron-hole pairs 


(rry) = Y AN dari Ger (Te) (16) 


kvc 


is a function of six coordinates, where r, and r, are the spatial 
coordinates of electron and hole. ø(r) are the Kohm-Sham or- 
bitals. We can define useful magnitudes from the exciton wave 
function. The weight of a transition v c is defined as the sum 


over all k 
Dy An: 


Analogously, we define the weight of each k, by summing 


over all transitions: 
wi = E 2, au 


The amplitude of electron-hole pairs that compose each ex- 
citon, as a function of the transition energy is 


(17) 
(18) 


Fw) = Y ale) so-oa. a9 
vck 


Finally, the optical absorption spectrum is the imaginary part 
of the dielectric function, e(hw), written as 


2:45 aleia | 


vk . L. écQ* = 


ex(hw) ec X, 


X 


ħw-T), (20) 


kvc 


where (óa|pilów) are the dipole matrix elements of the transi- 
tions v c. The vector p represents the light polarization. We 
restrict to light linearly polarized in the basal plane. The po- 
larization perpendicular to the basal plane of MoS% has a neg- 
ligible contribution to absorption for energies close to the band 
gap. 

Realistic results are only possible by performing an adequate 
convergence. In BSE calculations, we have to check carefully 
the number of valence and conduction bands, as well the k- 
point mesh. Coarse k-grids tend to overestimate the exciton 
binding energy. The building blocks of the BSE kernel, Ke}, are 
the screened Coulomb and the exchange interaction. Therefore, 
the dielectric function which enters in the Coulomb interaction 
has also to be converged with the number of bands and the k- 
point grid (see Refs. for details). 

The supercell geometry is also an important factor in BSE 
calculations. The long-range Coulomb interaction between 
replicas decreases slowly with distance. Consequently, GW and 
BSE corrections converge also slowly with the separation be- 
tween replicas [182]. In general, these corrections have 
opposite sign and partially cancel each other, and the total cor- 
rection of the band gap is close to the experimental value. How- 
ever, exact determination of the exciton binding energy requires 
overcoming this problem. An efficient technique is truncating 
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Figure 18: (a) Optical absorption with/without (solid/dashed lines) electron- 
hole interaction, together with the amplitude g(w) of excitons X4 and Xg. Ab- 
sorption thresholds are indicated by vertical solid/dashed lines. (b) Weigth "x 
of the exciton X4. 


the Coulomb potential (or Coulomb cut-off), simulating an in- 
finite distance between replicas [183]. In single-layer MoS», 
both GW band gap and exciton binding energy increases no- 
tably altough not their difference. A drawback of the Coulomb 
cut-off technique is the slower convergence with respect to the 
k-point grid [161]. 

Figure[18[a) shows a typical BSE calculation for single-layer 
MoS». Absoprtion is depicted with and without excitonic ef- 
fects (solid and dashed black curves). Electron-hole interaction 
down-shifts the absorption threshold and increases the absorp- 
tion coefficient. Amplitude functions g*^ (w) and g*#(w) show 
the typical profile of an exciton built from a transition between 
nearly parabolic bands. Contribution decays for increasing en- 
ergy, with the maximum close to the band edges. Panel (b) of 
Figure |18 represents the weight wie of the exciton A (for exci- 
ton B we obtain a similar result). The contribution is localized 
at K and the k-grid must be fine enough to describe accurately 
excitons A and B (145][156]. The distribution of the weight w^ 
reflects the importance of a proper convergence of the k-grid to 
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obtain realistic calculations. 

In the following, we discuss the calculations of Ref. [145]. 
The k-grids are 51 x 51 x 1 (for single- and double-layer), and 
21 x 21 x 3 for bulk. We have included the bands in the energy 
window from -3 to 3 eV. 

Figure shows the optical absorption for single-layer, 
double-layer and bulk MoS, (solid lines). For comparison, 
we have included the independent-particle absorption spectra, 
without the electron-hole interaction, obtained in the random- 
phase approximation (RPA, dashed lines). The starting point 
of BSE is the GW eigenvalues and the LDA wave functions. 
We have drawn the experimental optical absorption (dots). We 
have rigidly shifted the theoretical spectra to match with the 
experimental points. The discrepancy is around 0.2 eV, within 
the error margin of GW and BSE calculations. Nonetheless, 
BSE describes well the main trends of the spectra. We remark 
the agreement in the absorption threshold, where BSE spectra 
reproduce accurately the spin-orbit splitting. The theoretical 
absorption at high energies also matches with the experiments. 
These high-energy peaks come from transitions located around 
the I point. 

We can also make a comparative analysis of single-layer, 
double-layer and bulk MoS;. The spectra have the same line- 
shape at the absorption threshold, two peaks that correspond 
to excitons X4 and Xs, followed by a plateau. The differences 
arise from the exciton binding energy, which decreases with the 
number of layers. The reason of such decreasing is related to a 
larger dielectric screening in double-layer and bulk. The high- 
energy exciton exhibits a sharp peak in the case of single-layer 
MoS;, and it becomes difficult to distinguish in double-layer 
and bulk. Experimental observation confirms the latter result, in 
which we observe a broad absorption for bulk MoS;, in contrast 
to the relatively narrow peak of single-layer MoS;. Putting to- 
gether the theoretical and experimental data we can deduce that 
the interlayer interaction affects mainly exciton X4 and Xp. This 
conclusion agrees with the study of MoS2/WSz heterostructures 
of Ref. [173]: the optical spectra of MoS2/WS, is the addition 
of the spectra of independent single-layers of MoS. and WS», 
and not the combination of spectra modified by a strong inter- 
layer coupling. Inspecting the exciton wave functions we can 
get a better insight into the interlayer coupling. 

The intensity of the excitonic peaks is directly related to the 
spatial localization of the wave function. Figure [20] shows the 
exciton wave functions for the excitons (a) X4 and (b) Xy for 
the case of single-layer MoS». We have placed the hole at the 
Mo atom. This makes sense as the valence band states at K 
are composed primarily by Mo d-orbitals (see Table|5). The 
exciton X4 is extended more than 50 Á, in consonance with 
the localization at the momentum space (see Fig. [I8). Exciton 
Xp has an identical wave function (not shown here). On the 
contrary, the high-energy exciton is localized in a few unit cells, 
in less that 30 A. The exciton X y is built from a contribution of 
a larger set of k points than in the case of X4. We will discuss 
the features of this exciton below. 

By looking at the exciton wave function in a plane parallel to 
the vertical axis we can learn more details about the interlayer 
coupling. Figures poro and (d) show the lateral view of the 
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Figure 19: Optical spectra for single-, double-layer and bulk MoS2, BSE (solid 
lines) and RPA (dashed lines). The experimental data has been collected from 
several publications, red squares(43), red and blue circles[10], green circles(11) 
and green triangles [110]. 


X4 exciton in single-layer and bulk, respectively. The exciton 
density lies mainly on the Mo atoms for the single-layer and 
only a tiny fraction is outside the layer. Bulk MoS, does not 
exhibit big differences with respect to the single-layer. There- 
fore, excitons X4 and Xg remain in one layer, without coupling 
between layers, and optical transitions take place at the same 
layer. In other layered materials like hexagonal boron nitride, 
we find that excitons can spread over several layers [184]. 
Another distinctive optical property of single-layer MoS, is 
the exciton Xy, visualized by a sharp peak at high energy (2.75 
eV). Is this peak built from a single exciton or do we see the 
collective effect of many excitons? Figure 1](a) shows the ex- 
citonic peak Xy, together with each contribution (vertical lines). 
On the right side we have drawn the profile of the correspond- 
ing exciton wave function. Figure|21|(b) depicts the weigth wie 
of the first vertical line (in red). The characteristics of this exci- 
ton are radically different from the case of X4 and Xg excitons. 
First, the exciton is localized around T in the k space, forming a 
kind of hexagonal wheel. Second, defining the binding energy 
is ambiguous. We know that for bound excitons, like X4 and 
X5, the binding energy is defined as the difference between ex- 
citon energy and the band gap energy. The transition energies 
of the Xy exciton fall within the continuous of states, making 
difficult to define such binding energy. Third, the sharp peak is 
the collective contribution of several excitons with similar en- 
ergy and wave function, as we can see from the wave functions 
of Fig. [21] (b). The parallel transitions lead to a singularity in 
the density of states, and often the term Van-Hove exciton is 


used to denominate such peak [185] : 
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Figure 21: (a): Bethe-Salpeter spectra of MoS» single-layer together with the 
side view of the exciton wave functions (marked with a vertical line). (b) Band 
structure of single-layer MoS close to T and wave function of Xy exciton 
represented in k-space. 


Recently, several experiments have found fingerprints of the 
Van-Hove exciton, e.g., two-photon spectroscopy [43], pho- 
tocurrent spectroscopy [188], and light scattering spectra [189]. 
Some properties of the Van Hove exciton, under discussion 
nowadays, are the large electric field required to dissociate the 
exciton, and the spontaneous decay of Van Hove excitons into a 
free electron-hole pair. 

Additionally, single-layer MoS, exhibits a Rydberg-like ex- 
citon series in the optical spectra [156]. To capture these excited 
states the convergence of the k-grid requires up to 72x 72x 1. In 
comparison with the Rydberg series for a 2D hydrogen model, 
the excitation spectrum of single-layer MoS, is completely dif- 
ferent (see in supplementary material of Ref. (156]). The rea- 
son for this difference is the spatial variation of the dielectric 
function in MoS». The Rydberg-like series and its particular 
behaviour has been also observed in single-layer WS» [190]. 


5.2. Optical spectra of other transition metal dichalcogenides 


The rest of semiconductor transition metal dichalcogenides, 
MoSe2, WS», and WSe;, shares the interesting optical prop- 
erties of MoS2. They have similar lattice parameters, band gap 


and spin-orbit splitting [191] . However, band 
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Figure 20: Top view of excitons (a) X4 and (b) Xy in single-layer MoS». Lateral view of the exciton XA for (c) single-layer and (d) bulk MoS». 
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WS; WSe5 MoSe; 
a (Ang.) 3.127 3.250 3253 
E, (eV) 1.739 1.458 1.516 
Aso (meV) 406.8 457.2 185.6 


Table 10: Lattice parameters, band gap at K (£,), and spin-orbit splitting at K 
(Aso) of TMDs, obtained with DFT-LDA. 


gaps are different enough to generate band mismatch in het- 
erostructures. We can combine TMDs in vertical heterostruc- 
tures for quantum well growth [49] [173][194], with the purpose 
of selective confinement of photogenerated excitons. 


Figure shows the BSE spectra of single-layer MoSe», 
WS», and WSe» (solid lines) and the RPA spectra (dashed 
lines). We have included experimental data for WS2, and WSe;. 
On the right side we present a lateral view of the exciton wave 
functions. Starting points are the LDA calculations, including 
spin-orbit. We have used a 51 x 51 x 1 k-grid and we have 
included four valence and conduction bands. The static dielec- 
tric function is obtained with 60 bands. More accurate spectra 
require using self-consistent GW quasiparticle eigenvalues, as 
shown in Section [4] For introductory purposes, using LDA as 
starting point allows to describe the main physics of the optical 
properties of TMDs. 


All the spectra exhibit two well differentiated excitons, A and 
B, which come from transitions centered at the K point, anal- 
ogously to single-layer MoS2. The spin-orbit splitting deter- 
mines the separation between the peaks A and B. The theoreti- 
cal splitting agrees very well with the experiments. Exciton A is 
uniquely composed of the top of the valence band and the bot- 
tom of the conduction band. Accordingly, exciton B is mostly 
composed by the second valence band (K,2) and the conduction 
band with opposite spin. In comparison with the A exciton of 
single-layer MoS, we observe a similar spreading of the wave 
functions for the other TMDs. Evidently, the spin-orbit split- 
ting is much higher in compounds which include tungsten. In 
valley-physics, this has important consequences. The tuning of 
the excitation energy is crucial to obtain an excitonic population 
with certain polarization. Tuning in WS? and WSe» will in prin- 
ciple be easier due to the energy separation and the generation 
of a valley polarization will be more efficient. We will comment 
on this again in subsection[5.5] devoted to valley physics. 


The high-energy excitons (from H; to H3) show a strong spa- 
tial confinement, analogously to MoS2. However, they split in 
several and well differentiated peaks. Spin-orbit also splits the 
bands around I, and this results in the splitting of the excitons. 
The experimental spectra of WSe» agree with the calculation 
in the relative separation between the H; and H3 peaks. This 
latter compound presents the strongest spin-orbit coupling ef- 
fects, either for the bound excitons A and B, or the Van-Hove 
exctions Hs. In summary, selenium-based TMDs have an H- 
peak separation close to 0.5 eV and sulphur-based TMDs have 
a difference below 0.2 eV. 
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Figure 22: From top to bottom, optical spectra of single-layer MoSe2, WS», 
WSe», with and without excitonic effects (solid and dashed black line). The ex- 
citon wave functions on the right are associated with the vertical lines marked in 
the spectra. We have placed the hole 0.5 Bohr on top of Mo (W). Experimental 
data are extracted from Ref. [195]. 


5.3. Empirical models for excitons in MoS 


Analytical approaches are very useful in the study of excitons 
in single-layer MoS; and they can be easily expanded to explore 
excitonic complexes such as trions and biexcitons [40]. These 
approaches are also suitable to obtain the Rydberg excitonic 
series [190]. 

The tight-binding ansatz of the electronic wave function 
takes the orbital composition of the valence and conduction 
band states obtained from DFT results. Using the density ma- 
trix formalism, one obtains the analytical solution of the band 
structure close to the points K and K’ [174]. 


1 
Ge = +54 (Age)? + Alp fA), 


where € is the valley and A = vy, vj, c1, c, denotes band and 
spin. Momentum dependence is given by 


f(K) = 3 + 2cos(k,) + 4 cos(k,/2) cos( V3k,/2). 


The Taylor expansion simplifies the eigenvalue momentum 
dependence to a parabolic band structure 


and the solution of the model Hamiltonian gives the eigen- 
vectors, from which one can obtain the carrier-light matrix el- 
ements. Transitions at K and K’ are between Mo-d orbitals 
of the valence band and S-p orbitals of the conduction band. 
These transitions can be optically excited by circularly polar- 
ized light. The right-handed circularly polarized light will ex- 
cite states at K and the left-handed light at K’, allowing a valley 
selection. This is the cornerstone of valley physics, which will 
be presented in subsection [5.5] Figure 23|shows the results of 
Ref. for the matrix elements and the absorption spec- 
tra for negative (a,c) and positive (b,d) light polarization. Fig. 
IQ shows our own calculations of the optical matrix elements 
for linear polarization. In the latter case we excite the valleys 
K and K' with the same probability. Agreement between DFT 
and the model Hamiltonian is excellent and justifies the use of 
the analytical approach. 

Excitonic effects are included introducing the Coulomb inter- 
action into the model Hamiltonian [174]. In the framework of 
the semiconductor Bloch equation, from microscopic polariza- 
tion, one can obtain an analytical expression of the absorption 
coefficient (for a complete derivation, see Ref. (196]). In addi- 
tion, the effects of the substrate on the exciton binding energy 
can be quantified by a proper choice of the dielectric constant. 

The application of the model shows a binding energy of the 
X4 exciton for free-standing single-layer MoS, of 860 meV. 
The binding energy decreases to 455 meV on top of a silicon 
oxide substrate. Among limitations of this analytical approch 
are the calculation of the relative intensity of the X4 and Xg 
peaks, which can be attributed to higher-order effects beyond 
the Hartree-Fock approximation. The prediciton of the high- 
energy excitons will also require a much more complicated re- 
formulation of the model. 
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Figure 23: (a) and (b): Optical matrix elements for negative (o—) and positive 
(a+) circular polarization. The corresponding spectra are represented in (c) 
and (d) (reprinted with permission from Ref. [174]. Copyright (2014) by the 
American Physical Society). (e) Optical matrix elements calculated with Yambo 
for linear light polarization. 


5.4. Experimental determination of the band gap and exciton 
binding energy in MoSe» 


The experimental determination of electronic band gap and 
exciton binding energy requires at least two techniques for de- 
termining univoquely these magnitudes. We have seen that 
the electronic band gap is related to the single-particle exci- 
tation. In addition, the binding energy is the difference be- 
tween the electronic and the optical band gap. In Ref. [197], 
Ugeda et. al. have used a combination of experiments and 
theory to give an accurate value for the binding energy of 
single-layer MoSe;. The measurements by scanning tunnel- 
ing spectroscopy (STS) have measured the electronic gap, and 
photoluminescence (PL) has defined the optical gap. The re- 
ported values are E, = 2.18 + 0.04 eV (electronic band gap) 
and Ej, = 1.63 + 0.01 eV, what gives a binding energy of 
0.55 + 0.04. Experimental findings are well supported by GW 
and BSE calculations, taking into account the incidence of the 
substrate (bilayer graphene). As mentioned, the substrate in- 
creases the dielectric constant and reduces the electronic and 
optical band gaps. However, a full GW and BSE calculation 
of the system MoSe; plus substrate would be computationally 
very complex. Alternatively, authors have made apart the calcu- 
lation of the substrate dielectric screening. The MoSe; contri- 
bution is obtained in the random phase approximation, includ- 
ing local fields. Afterwards, the MoSe» and substrate contri- 
butions are merged in the Bethe-Salpeter equation. The bind- 
ing energy reduces from a value of 0.65 eV to 0.52 eV, with 
an uncertainty of 0.10 eV, in fair agreement with experimental 
values. Therefore, a wise treatment of substrate influence on 
electronic and optical properties appears as an important aspect 
in calculations aiming of having predictive character. 


5.5. Spin-orbit interaction and valley physics in MoS 


The lack of inversion symmetry and the strong spin-orbit in- 
teraction in single-layer MoS» lead to what is called valley- 
physics. The valley index denotes the momentum of the valence 
band state, K or K', and the spin (up or down). As shown above, 
excitons A and B can be generated exclusively from the valley 
K or K’ by selecting the appropiate light polarization. Among 
the potential uses and research related with this new concept 
are information transport by means of a new carrier, defined 
in terms of the valley index, or the generation of a valley-Hall 
effect [38] [36] [37] [198]. 

In the case of single-layer MoS2 , spin-orbit coupling splits 
the valence band maximum at K and K’ by ~ 150 meV. More- 
over, the point group symmetry Ds; does not have inversion 
symmetry. Under these conditions, Kramer’s degeneracy states 
that E;(k) = E,(—k). In the case of the valence band states at K, 
we find the relation E;(K) = E,(K’). The spin of the valence 
band states K,; and Ky» flips under change of inequivalent K 
points. In other words, valence band edge valley and spin are 
coupled. Figure Baa) illustrates the spin composition of the 
valence band edges at K and K’. The spin-orbit splitting is neg- 
ligible (a few meV) for conduction states and we can consider a 
two-fold degenerate state. Following Ref. [35], we can assign 
the following total angular momenta, mj, in the z-direction to 
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Figure 24: (a) Illustration of the valence and conduction bands at the mathb f K 
points of the Brillouin zone, including the spin-orbit interaction (reprinted with 
permission of Ref. [35], copyright (2014) by the American Physical Society). 
(b) (reprinted by permission from Macmillan Publishers Ltd: Nature Nanotech- 
nology, Ref. , copyright 2012. 


the conduction and valence band states. The conduction band 
state at K is energy degenerate with m; = —3/2 and m; = —1/2. 
In the case of valence band states, K,; and K,? have m; = —1/2 
and m; = 1/2 respectively, split in energy by ~ 150 meV. The 
values for the total angular momentum of the valence and con- 
duction states at K’ are obtained by multiplying the ones at K 
by -1. 

Exposing the crystal to circularly polarized light of mo- 
mentum +1 will promote a valence electron with momentum 
m; = +1/2 to the conduction state of m; = +3/2, or from the 
valence state m; = —1/2 to the conduction state m; = +1/2, 
both transitions taking place in the valley K'. These transitions 
are consequently forbidden at valley K. By tuning the excita- 
tion energy one can precisely select which couple of valence- 
conduction band states are participating and from which valley. 
Thus, light with the energy of the band gap and momentum +1 
will only promote valence electrons with m; = +1/2. In this 
way, one can create a stable population of electron-hole pairs 
with defined valley index. In other words, a “valley polariza- 
tion" is generated. This has been demonstrated theoretically 
and experimentally in Refs. (38][36][35]. The persistence of this 
valley polarization is related to the conservation of spin. Hence, 
the flipping of valley index implies a spin flip, which is only 
possible by scattering with magnetic impurities or by means of 
relaxation via intra- and inter-valley scattering [35][199]. 


Nowadays, the selective generation of valley polarizations is 
presented as a way to develop the so-called valleytronics and the 


stability of the valley polarization is being investigated. Tech- 
niques such as ultra-fast spectroscopy make it possible to trace 
the time evolution of valley polarization. Although theory pre- 
dicts a long lifetime for the valley polarization, recent stud- 
ies observe a non-negligible decoherence [200], which points 
towards other reasons than magnetic impurities as the origin 
of the spin-flip. Experimental studies have been undertaken 
to identify the reasons of decoherence of valley polarizations 
[199]. The analysis of the time-dependent 
optical response of single-layer MoS; supported in several sub- 
strates or suspended in vacuum can give a hint of how 
the environment can alter the electronic structure and hence the 
valley polarization. Other causes of decoherence can be the in- 
terplay between intra-valley and inter-valley scattering. Here 
the detuning of the valley polarization can take place via inter- 
mediate transitions through the T point [199]. The proliferation 
of experimental results using time-dependent spectroscopy is 
increasing our knowledge of the electronic and optical proper- 
ties of MoS, and other two-dimensional materials. 
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6. Conclusions 


This review has summarized the theoretical and computa- 
tional description of the fundamental spectroscopic properties 
of MoS, in its single-layer, few-layer, and bulk form. We have 
focussed on MoS, but many of the findings are similarly valid 
for the other semiconducting transition-metal dichalcogenides 
(TMDs) of 2H polytype. We have summarized the numerous 
published investigations that report partly controversial find- 
ings. We have analysed the possible reasons for the discrep- 
ancies by performing comprehensive density functional theory 
calculations on different levels of approximation. For the ge- 
ometry, we have investigated the inclusion of Van der Waals 
interactions. For the quasiparticle gap, we have compared the 
results obtained by DFT with different hybrid functionals with 
the result obtained by many-body perturbation theory. Exci- 
tonic effects in the optical properties are described with the 
Bethe-Salpeter approach. Thereby, our review provides a gen- 
eral idea on the most important computational issues that arise 
when studying TMDs. We hope that it gives stimulation to the 
scientific community to achieve accurate and converged results. 

The structural and vibrational properties are well described 
by density functional theory approaches including Van der 
Waals interactions. Excellent results are obtained for the lattice 
parameters and the calculated phonon frequencies that agree 
well with different experimental data from Raman spectroscopy 
and neutron scattering. In this context, the open question about 
the anomalous Davydov splitting has been explained in terms 
of many neighbors interaction between Mo-S atoms of differ- 
ent layers. The anomalous trend of the E», mode as a function 
of layer number is a consequence of the renormalization of the 
atomic distances due to the free surface[62]. 

There is a variety of results in the literature for the band struc- 
tures and band gaps based on different levels of computational 
approach. The spread of results is connected to the inherent 
problem of local and semilocal exchange correlation function- 
als to yield accurate band gaps and can be overcome by includ- 
ing nonlocal exchange (hybrid functionals) or using the GW ap- 
proximation. The latter is commonly applied in a one-shot man- 
ner (GoWo) on top of DFT wave functions. Different schemes 
of self-consistency lead to a large spread in the quasiparticle 
gaps. Reliable band structures and relative positions between 
the two lowest conduction band extrema K, and T,, but slightly 
overestimated band gaps, are obtained by self-consistent quasi- 
particle GW calculations provided that one starts from a fully 
optimized crystal structure. Particularly important in the case 
of single-layer MoS, is the convergence of the calculations with 
respect to the k-point grid, the number of virtual states, and the 
vacuum layer in the slab approach. 

As for most other layerd materials, excitonic effects are very 
pronounced in the optical properties of MoS2. The absorption 
spectra of mono-, bilayer, and bulk MoS, display three pro- 
nounced peaks. The excitons are calculated with the Bethe- 
Salpeter equation based on the full spinorial wave functions 
in order to include the effects of spin-orbit interaction. In all 
single-layer, double-layer and bulk MoS;, there is a pronounced 
splitting between A and B excitons which can be traced back 
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to the splitting of the valence band maximum at the high- 
symmetry point K. For single-layer MoS;, the splitting is due 
to the strong spin-orbit splitting of the valence band maximum 
at K, for for double-layer and bulk MoS,, the splitting is due to 
the interlayer interaction. Interestingly, the brightest exciton of 
single-layer MoS, is not found at the absorption threshold, but 
at higher energies (around 3 eV). This exciton, also called “Van 
Hove exciton" stems from a large joint density of states due to 
parallel conduction and valence bands around I. It plays an 
important role in the resonant Raman spectroscopy of various 
semiconducing transition-metal dichalcogenides. 
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